{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* By: Illya Barziy\n",
    "* Email: illyabarziy@gmail.com\n",
    "* Reference: __A Robust Estimator of the Efficient Frontier__ _by_ Marcos Lopez de Prado"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The Nested Clustered Optimization algorithm (NCO)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This description of the algorithm and the realizations in the mlfinlab package are based on the paper by _Marcos Lopez de Prado_ __A Robust Estimator of the Efficient Frontier__  [available here](https://papers.ssrn.com/abstract_id=3469961)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Introduction"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Convex optimization solutions in the problems of portfolio optimization tend to be unstable. The instability comes from two sources:\n",
    "- Noise in the input variables\n",
    "- The signal structure that magnifies the estimation errors in the input variables\n",
    "\n",
    "The NCO algorithm tackles both sources of instability."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Problem Statement"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This part is a direct quote from the paper __A Robust Estimator of the Efficient Frontier__ with minor comments.\n",
    "\n",
    "Say, we have a system with $N$ random variables, the expected value of draws of the variables are $\\mu$ and the covariance matrix is $V$. We want to compute the vector $w$ that minimizes the variance of the system, where the variance is $w'Vw$, subject to $w'a$. Here $a$ characterizes the optimal solution. \n",
    "\n",
    "So, the problem is:\n",
    "\n",
    "$$min_{w}\\frac{1}{2}w'Vw$$\n",
    "\n",
    "$$s.t.: w'a = 1 $$\n",
    "\n",
    "In financial application, the optimal allocation $w^*$ known as the maximum Sharpe ratio, when $a = \\mu$ (as it maximizes $\\frac{w'\\mu}{\\sqrt{w'Vw}}$).\n",
    "\n",
    "The optimal allocation $w^*$ known as the minimum variance portfolio, when $a = 1_{N}$, where $1_{N}$ is the vector of ones of size $N$ (as it minimizes the $w'Vw$ - variance).\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Using the convex optimization methods, the optimal solution to the problem is\n",
    "\n",
    "$$w^* = \\frac{V^{-1}a}{a'V^{-1}a}$$\n",
    "\n",
    "This is the Convex Optimization Solution (CVO). It allows finding the optimal weights $w$, when the real $V$ and $a$ of the variables are known."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As the input variables $V$ and $a$ are typically unknown, the estimations are used, which leads to unstable solutions where a small change of inputs will cause bix changes of $w^*$.\n",
    "\n",
    "The estimators are obtained from observations (draws) of $N$ variables, or instruments in a portfolio."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Instability caused by noise"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This part is a direct quote from the paper __A Robust Estimator of the Efficient Frontier__ with minor comments.\n",
    "\n",
    "We have a matrix of i.i.d. random observations $X$, of size $TxN$. The underlying process is generating the observations with zero mean and $\\sigma^2$ variance. \n",
    "\n",
    "The matrix $C = T^{-1}X'X$ has eigenvalues $\\lambda$ that asymptotically converge (as $N \\rightarrow +\\infty$ and $T \\rightarrow +\\infty$ with $1 < \\frac{T}{N} < +\\infty$) to the Marcenko-Pastur probability density function:\n",
    "\n",
    "$$ f[\\lambda] =\n",
    "  \\begin{cases}\n",
    "    \\frac{T}{N}\\frac{\\sqrt{(\\lambda_+-\\lambda)(\\lambda-\\lambda_-)}}{2\\pi\\lambda\\sigma^2} & \\quad \\text{if } i \\in [\\lambda_-,\\lambda_+] \\text{,}\\\\\n",
    "    0  & \\quad \\text{if } i \\notin [\\lambda_-,\\lambda_+] \\text{.}\n",
    "  \\end{cases}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Where the maximim expected eigenvalue is \n",
    "\n",
    "$$\\lambda_+ = \\sigma^2(1 + \\sqrt{\\frac{N}{T}})^2$$\n",
    "\n",
    "and the minimum expected eigenvalue is \n",
    "\n",
    "$$\\lambda_- = \\sigma^2(1 - \\sqrt{\\frac{N}{T}})^2$$\n",
    "\n",
    "When $\\sigma^2 = 1$, then $C$ is the correlation matrix associated with $X$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The Marcenko-Pastur distribution explains why empirical covariance matrices contain substantial amounts of noise. \n",
    "\n",
    "In many practical applications, $\\frac{N}{T}\\rightarrow 1$, thus the covariance's eigenvalues span a wide range $[\\lambda_-,\\lambda_+]$. Because $\\lambda_- \\rightarrow 0$ as $\\frac{N}{T}\\rightarrow 1$, the determinant of $V$ approaches zero and $V^{-1}$ cannot be estimated robustly, this is the source of unstable solutions $w^*$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The process of de-noising the covariance matrix that is used in the code examples in the end is described in a paper by _Potter M._, _J.P. Bouchaud_, _L. Laloux_ __“Financial applications of random matrix theory: Old laces and new pieces.”__  [available here](https://arxiv.org/abs/physics/0507111).\n",
    "\n",
    "It is briefly explained later in the notebook."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Instability caused by signal"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This part is a direct quote from the paper __A Robust Estimator of the Efficient Frontier__ with minor comments.\n",
    "\n",
    "Certain covariance structures can make the optimal solution unstable. Looking at the correlation matrix between two variables:\n",
    "\n",
    "$$C =\n",
    "  \\begin{bmatrix}\n",
    "    1 & \\rho\\\\\n",
    "    \\rho & 1 \n",
    "  \\end{bmatrix}$$\n",
    "\n",
    "Where the $\\rho$ is the correlation between their outcomes. Matrix $C$ can be diagonalized as $CW = W\\Lambda$, where"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\\Lambda = \n",
    "  \\begin{bmatrix}\n",
    "    1 & +\\rho \\\\\n",
    "    1 & -\\rho \n",
    "  \\end{bmatrix}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$W = \n",
    "  \\begin{bmatrix}\n",
    "    \\frac{1}{\\sqrt{2}} & \\frac{1}{\\sqrt{2}} \\\\\n",
    "    \\frac{1}{\\sqrt{2}} & -\\frac{1}{\\sqrt{2}} \n",
    "  \\end{bmatrix}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The inverse of $C$ is\n",
    "\n",
    "$$C^{-1}=W\\Lambda^{-1}W'=\\frac{1}{|C|} =\n",
    "  \\begin{bmatrix}\n",
    "    1 & -\\rho \\\\\n",
    "    -\\rho & 1 \n",
    "  \\end{bmatrix}$$\n",
    "\n",
    "Where $|C|$ is the determinant of $C$, $|C|=(1+\\rho)(1-\\rho)=1-\\rho^2$.\n",
    "\n",
    "We can see that more $\\rho$ derives from zero, the bigger one eigenvalue becomes relative to the other, causing the determinant of $C$ to approach zero, which makes the values of $C^{-1}$ explode."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "When the correlation matrix is an identity matrix, the eigenvalue function is a horizontal line. Outside that ideal case, at least one subset of variables exhibits greater correlation among themselves than to the rest, forming a cluster within the correlation matrix. \n",
    "\n",
    "When $K$ variables form a cluster, they are more heavily exposed to a common eigenvector, which implies that the associated eigenvalue explains a greater amount of variance. \n",
    "\n",
    "But because the trace of the correlation matrix is exactly $N$, that means that an eigenvalue can only increase at the expense of the other $K-1$ eigenvalues in that cluster, resulting in a condition number greater than 1. Consequently, the greater the intra-cluster correlation is, the higher the condition number becomes. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### De-noising function"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The function provided below for de-noising the covariance works as follows:\n",
    "- The given covariance matrix is transformed to the correlation matrix.\n",
    "- The eigenvalues and eigenvectors of the correlation matrix are calculated.\n",
    "- Using the Kernel Density Estimate algorithm a kernel of the eigenvalues is estimated.\n",
    "- The Marcenko-Pastur pdf is fitted to the KDE estimate using the variance as the parameter for the optimization.\n",
    "- From the obtained Marcenko-Pastur distribution, the maximum theoretical eigenvalue is calculated using the formula from the \"Instability caused by noise\" part.\n",
    "- The eigenvalues in the set that are above the theoretical value are all set to their average value. For example, we have a set of 5 sorted eigenvalues ($\\lambda_1$...$\\lambda_5$), 2 of which are above the maximum theoretical value, then we set $\\lambda_4^{NEW} = \\lambda_5^{NEW} = \\frac{\\lambda_4^{OLD} + \\lambda_5^{OLD}}{2}$\n",
    "- The new set of eigenvalues with the set of eigenvectors is used to obtain the new de-noised correlation matrix.\n",
    "- The new correlation matrix is then transformed back to the new de-noised covariance matrix."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is how the de-noising function from mlfinlab package can be used:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import mlfinlab as ml\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The de-noised covariance matrix is:\n",
      "[[ 0.01        0.00267029 -0.00133514]\n",
      " [ 0.00267029  0.04       -0.00438387]\n",
      " [-0.00133514 -0.00438387  0.01      ]]\n"
     ]
    }
   ],
   "source": [
    "# A class that has the de-noising function\n",
    "risk_estimators = ml.portfolio_optimization.RiskEstimators()\n",
    "\n",
    "# Covariance matrix to de-noise\n",
    "cov_matrix = np.array([[0.01, 0.002, -0.001],\n",
    "                       [0.002, 0.04, -0.006],\n",
    "                       [-0.001, -0.006, 0.01]])\n",
    "\n",
    "# Relation of number of observations T to the number of variables N (T/N)\n",
    "tn_relation = 50\n",
    "\n",
    "# The bandwidth of the KDE kernel\n",
    "kde_bwidth = 0.25\n",
    "\n",
    "# Finding the de-noised covariance matrix\n",
    "cov_matrix_denoised = risk_estimators.denoise_covariance(cov_matrix, tn_relation, kde_bwidth)\n",
    "\n",
    "# Outputting the result\n",
    "print('The de-noised covariance matrix is:')\n",
    "print(cov_matrix_denoised)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As we can see, the main diagonal hasn't changed, but the other covariances are different. This means that the algorithm has changed the eigenvalues of the correlation matrix."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Also functions to:\n",
    "- transform covariance matrix into correlation matrix \n",
    "- transform correlation matrix into covariance matrix \n",
    "\n",
    "are available in the mlfinlab package:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The correlation matrix is:\n",
      "[[ 1.   0.1 -0.1]\n",
      " [ 0.1  1.  -0.3]\n",
      " [-0.1 -0.3  1. ]]\n",
      "The covariance matrix calculated back is:\n",
      "[[ 0.01   0.002 -0.001]\n",
      " [ 0.002  0.04  -0.006]\n",
      " [-0.001 -0.006  0.01 ]]\n",
      "Exactly the same as the original one:\n",
      "[[ 0.01   0.002 -0.001]\n",
      " [ 0.002  0.04  -0.006]\n",
      " [-0.001 -0.006  0.01 ]]\n"
     ]
    }
   ],
   "source": [
    "# Transforming our covariance matrix to a correlation matrix\n",
    "corr_matrix = risk_estimators.cov_to_corr(cov_matrix)\n",
    "\n",
    "# Outputting the result\n",
    "print('The correlation matrix is:')\n",
    "print(corr_matrix)\n",
    "\n",
    "# The standard deviation to use when calculating the covaraince matrix back\n",
    "std = np.array([0.1,0.2,0.1])\n",
    "\n",
    "# And back to the covariance matrix\n",
    "cov_matrix_again = risk_estimators.corr_to_cov(corr_matrix, std)\n",
    "\n",
    "# Outputting the result\n",
    "print('The covariance matrix calculated back is:')\n",
    "print(cov_matrix_again)\n",
    "print('Exactly the same as the original one:')\n",
    "print(cov_matrix)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### NCO and CVO functions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The Convex Optimization Solution (CVO) solves the problem described in the \"Problem Statement\" part using the provided formula:\n",
    "\n",
    "$w^* = \\frac{V^{-1}a}{a'V^{-1}a}$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The NCO algorithm is following these steps:\n",
    "\n",
    "- Gets the covariance matrix of the outcomes as an input (and the vector of means if the target is to maximize the Sharpe ratio).\n",
    "- Transforms the covariance matrix to the correlation matrix and calculate the distance matrix based on it.\n",
    "- Clusters the covariance matrix into subsets of highly-correlated variables.\n",
    "- Computes the optimal weights allocation for every cluster (using the CVO).\n",
    "- Reduces the original covariance matrix to a reduced one - where each cluster is represented by a single variable.\n",
    "- Computes the optimal weights allocation for the reduced covariance matrix (using the CVO).\n",
    "- Computes the final allocations as a dot-product of the allocations between the clusters and inside the clusters.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is how the NCO and the CVO from mlfinlab package can be used:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The optimal weights with minimum variance using the NCO are:\n",
      "[[0.43298969]\n",
      " [0.13402062]\n",
      " [0.43298969]]\n"
     ]
    }
   ],
   "source": [
    "# A class that has the NCO and CVO functions\n",
    "nco = ml.portfolio_optimization.NCO()\n",
    "\n",
    "# Covariance matrix\n",
    "cov_matrix = np.array([[0.01, 0.002, -0.001],\n",
    "                       [0.002, 0.04, -0.006],\n",
    "                       [-0.001, -0.006, 0.01]])\n",
    "\n",
    "# Vector of ones (as our goal is minimizing the variance)\n",
    "mu_vec = np.array([1, 1, 1]).reshape(-1, 1)\n",
    "\n",
    "# Maximum number of clusters to use in the correlation matrix clustering\n",
    "max_num_clusters = 2\n",
    "\n",
    "# Finding the optimal weights using the NCO\n",
    "w_nco = nco.allocate_nco(cov_matrix, mu_vec, max_num_clusters=max_num_clusters)\n",
    "\n",
    "# Outputting the result\n",
    "print('The optimal weights with minimum variance using the NCO are:')\n",
    "print(w_nco)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now using the simple CVO algorithm:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The optimal weights with minimum variance using the CVO are:\n",
      "[[0.37686939]\n",
      " [0.14257228]\n",
      " [0.48055833]]\n"
     ]
    }
   ],
   "source": [
    "# Finding the optimal weights using the CVO\n",
    "w_cvo = nco.allocate_cvo(cov_matrix, mu_vec)\n",
    "\n",
    "# Outputting the result\n",
    "print('The optimal weights with minimum variance using the CVO are:')\n",
    "print(w_cvo)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The results are different, so the NCO has clustered the correlation matrix. We can see in later examples if this produces a more robust result."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### MCOS function"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The Monte Carlo Optimization Selection (MCOS) algorithm calculates the NCO allocations and a simple optimal allocation for multiple simulated pairs of mean vector and the covariance matrix to determine the most robust method for weights allocations for a given pair of means vector and a covariance vector. Basically, it will show if the NCO or the CVO method is more robust.\n",
    "\n",
    "However, the MCOS may support other optimization methods and compare their robustness.\n",
    "\n",
    "The steps of the MCOS algorithm are:\n",
    "- Gets the covariance matrix and the means vector of the true outcomes as an input (along with the simulation parameters to use).\n",
    "- Draws the empirical covariance matrix and the empirical means vector based on the true ones.\n",
    "- If the kde_bwidth parameter is given, the empirical covariance matrix is de-noised.\n",
    "- Based on the min_var_portf parameter, either the minimum variance or the maximum Sharpe ratio is targeted in weights allocation.\n",
    "- CVO is applied to the empirical data to obtain the weights allocation.\n",
    "- NCO is applied to the empirical data to obtain the weights allocation.\n",
    "- Based on the original covariance matrix and the means vector a true optimal allocation is calculated.\n",
    "- For each weights estimation in a method, a standard deviation between the true weights and the obtained weights from the simulation is calculated.\n",
    "- The error associated with each method is calculated as the mean of the standard deviation across all estimations for the method."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is how the MCOS from mlfinlab package can be used:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The optimal weights with minimum variance using the NCO from simulations are:\n",
      "          0         1         2         3\n",
      "0  0.181431  0.326042  0.165669  0.326857\n",
      "1  0.221613  0.394977  0.204978  0.178432\n",
      "The optimal weights with minimum variance using the CVO from simulations are:\n",
      "          0         1         2         3\n",
      "0  0.165431  0.344039  0.168434  0.322096\n",
      "1  0.202771  0.427418  0.212893  0.156919\n"
     ]
    }
   ],
   "source": [
    "# Covariance matrix\n",
    "cov_mat = np.array([[1, 0.1, 0.2, 0.3],\n",
    "                    [0.1, 1, 0.1, 0.2],\n",
    "                    [0.2, 0.1, 1, 0.1],\n",
    "                    [0.3, 0.2, 0.1, 1]])\n",
    "\n",
    "# Vector of means\n",
    "mu_vec = np.array([0, 0.1, 0.2, 0.3])\n",
    "\n",
    "# Number of observations to use to obtain empirical covariance matrix and the means vector\n",
    "num_obs = 100\n",
    "\n",
    "# Number of simulations to do in the MCOS\n",
    "num_sims = 2\n",
    "\n",
    "# The bandwidth of the KDE kernel.\n",
    "kde_bwidth = 0.25\n",
    "\n",
    "# Flag if the goal is the minimum variance\n",
    "min_var_portf = True\n",
    "\n",
    "# Flag if we want to use the Ledoit-Wolf shrinkage procedure\n",
    "lw_shrinkage = False\n",
    "\n",
    "# For the same result output\n",
    "np.random.seed(1)\n",
    "\n",
    "# Finding the optimal weights for minimum variance\n",
    "w_cvo, w_nco = nco.allocate_mcos(mu_vec, cov_mat, num_obs, num_sims, kde_bwidth, min_var_portf, lw_shrinkage)\n",
    "\n",
    "# Outputting the result\n",
    "print('The optimal weights with minimum variance using the NCO from simulations are:')\n",
    "print(w_nco)\n",
    "print('The optimal weights with minimum variance using the CVO from simulations are:')\n",
    "print(w_cvo)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And the outputs can be analyzed based on the mean of the standard deviation with the true weight allocations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The mean of the standard deviation with the true weight for NCO is:\n",
      "0.0524762395749268\n",
      "The mean of the standard deviation with the true weight for CVO is:\n",
      "0.05839900155706323\n"
     ]
    }
   ],
   "source": [
    "# Finding the errors in estimations\n",
    "err_cvo, err_nco = nco.estim_errors_mcos(w_cvo, w_nco, mu_vec, cov_mat, min_var_portf)\n",
    "\n",
    "# Outputting the result\n",
    "print('The mean of the standard deviation with the true weight for NCO is:')\n",
    "print(err_nco)\n",
    "print('The mean of the standard deviation with the true weight for CVO is:')\n",
    "print(err_cvo)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Even on this small example, we can see that the NCO has performed better than the CVO method (the error metric is smaller)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Sample Data Generating function"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can create a random vector of means and a random covariance matrix to properly test the algorithms described above. \n",
    "\n",
    "These created covariance matrix and vector of means have similar characteristics to what the real securities inside a portfolio may pose. \n",
    "\n",
    "The elements are divided into clusters. The elements in clusters have a given level of correlation. The correlation between the clusters is set at another level. This structure is created in order to test the NCO and MCOS algorithms."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is how the Sample Data Generating function from mlfinlab package can be used:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The generated vector of means is:\n",
      "[[-0.00089985]\n",
      " [ 0.2608328 ]\n",
      " [ 0.45695781]\n",
      " [ 0.0015997 ]]\n",
      "The generated covariance matrix is:\n",
      "       0      3      2      1\n",
      "0  0.090  0.000  0.000  0.027\n",
      "3  0.000  0.090  0.027  0.000\n",
      "2  0.000  0.027  0.090  0.000\n",
      "1  0.027  0.000  0.000  0.090\n"
     ]
    }
   ],
   "source": [
    "# Number of blocks to have in a modeled matrix\n",
    "num_blocks = 2\n",
    "\n",
    "# The size of each block in a modeled matrix\n",
    "block_size = 2\n",
    "\n",
    "# Correlation inside a block\n",
    "block_corr = 0.3\n",
    "\n",
    "# Correlation between the clusters\n",
    "std = 0.3\n",
    "\n",
    "# Finding the random vector of means and covariance matrix\n",
    "mu_vec, cov_matrix = nco.form_true_matrix(num_blocks, block_size, block_corr, std)\n",
    "\n",
    "# Outputting the result\n",
    "print('The generated vector of means is:')\n",
    "print(mu_vec)\n",
    "print('The generated covariance matrix is:')\n",
    "print(cov_matrix)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This covariance vector and matrix of means can be used to test the NCO and the CVO methods."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Output visualization"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's use all of the described above functions on a bigger set of generated data and visualize the results."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "import seaborn as sns\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.axes._subplots.AxesSubplot at 0x20f69b90ec8>"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWQAAAD9CAYAAACLBQ0fAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3df3wU9bn//de1mxDAYAwhIZAgUHEXsT21oqktFRVvFSlCexctglZFmtYj941VWgv65Vi09uuv0+oNnoKUKphWOfZ7KtqgtlIsCipIVX4ZC4KyECSAIUT5kcxe9x+7xE0I2Qls3JnlevqYhzszn5l9Z4FrP/nMZ2dFVTHGGJN+gXQHMMYYE2MF2RhjPMIKsjHGeIQVZGOM8QgryMYY4xFWkI0xxiM8WZDD4fDwcDhcFQ6HN4bD4Z+3sn9oOBxeHQ6HG8Ph8JgW+04Nh8MvhcPhDeFweH04HO7XgVGHA1XARuCInEAO8HR8/xvA4SydgN8Da4B3gAsTjukEzAHeB94DvncC522vZD+fF/kts9/y+ornCnI4HA4Cs4DLgUHA1eFweFCLZh8B1wN/aOUU84EHqqqqzgDKgJ0dFPWInPH/J7oR+AQYAPwauC++/Yfx/38FuAR4iM//LO6IZw7Fz/fKCZq3vdz8fF7jt8x+y3vcRGS4iFSJyEYROeINSESGishqEWkUkTEt9p0qIi+JyAYRWS8i/ZI9X9KCLCJlInJu/PEgEblVREa4/5HarQzYWFVV9UFVVdUh4ClgdGKDqqqqLVVVVe8C0cTt8cKdVVVV9dd4u/qqqqrPOjIn8AHQas74+hPxx88AFwNC7C/xy/HtO4Fa4Jz4+gTgV/HHUWDXCZq3vdz8fF7jt8x+y3tcROSINyARaXfnUFVddw7bLMgi8h/AI8B/icivgJlALvBzEbkj2cmPUQmwNWE9Et/mRgioDYfD/yccDv8zHA4/EO9xdwQ3ORPbNAJ7gQJiv/aPBrKA/sBgoA9wSrzt3cBq4L+Bnido3vY6nr836eK3zH7Le7zKgI2q+oGqtvoGpKpbVPWIzmG8cGep6l/j7epVNWnnUNr66LSIrAHOIja2uAMoVdU6EekCvKGq/3aU48qBcoBHH7pn8MQfXJ0sR5MXlyzjtTfeYsbUWwBY9MLLrF1fxbRb//2Itnfc8xAXDCnj0ovOB+Clvy9j+q9+w3//fia9ehYxZfqvOP8b5/K9Ky5z/fx7xkxw1S7nggvIKSuj7oEHAOh8ySVkn3EG+x55pKlNwe9/zyc/+xnRmprYekUFe266Cf30U3J//GM6fe1rODt2IFlZ7H/uOQ6tXUvRs89SO306B//xD7peeSVZp59O3b33us7vl7wly/913D9Tou99bySXXnIBP/rxTwEYP/57nHvOWdzyk/+V0udJJb9l/iLyNh7aJsd7joZdH7i+H0SnwtN+RLxWxc1R1TkA8SGI4ao6Mb5+LfB1VZ3U8jwi8jjwvKo+E1//DjCR2G8S/YG/AT9XVaetPFlJ8jbGT/CZiGxS1ToAVd0vItGjHRT/geZA+14cgJ5FPdixs6Zp/eOduyjsUeDu2MIeDAydRp+SXgAMG/oN3l33HuC+ILsVrakhUFjYtB4oLMTZ1fy3daemhmBhYazABYMEcnPRujoA6mfNamqXP3MmjZEIuncvun8/B5ctA+DA0qXkj0jN6JDf8rbXtkg1fUp7N62XlvSiuvrjtGRxy2+Z/ZbXjcRa1YrW3hzc1rMs4Hzga8SGNZ4mNrTxu7YOSjaGfEhEusYfD25KKZJHiy56qnx5YIiPItuJbN9BQ0MDi19+hYu+dZ67Y88IUbevnj2f1ALw5lvvcFq/UzsiJg1VVQRLSwkUF0NWFp2HDePg8uXN2hxcvpzOw4cDsR7qodWrYztycqBzZwA6DR4MjoPz4YexY1asIPuss5r2Nca3n2h522vlqrcZMKA//fr1ITs7m6uuGs1zz7+Ulixu+S2zb/JGHfdL2yLEhuYOKwW2u0wRAf4ZH+5oBP4MnJ3soGQ95KGqehBAVRMLcDZwnctg7ZKVFWTaT27iR7feieM4fHfkpQz4Ul9mPjafMweGuOj881izoYpbpt5N3b56lr72BrPmPsmzFbMJBoNMuXkiN06eCgqDwgMYM2p4R8QEx2Hfww+T/8ADEAhwYPFinC1bOOmGG2isquLg8uXsr6wkb9o0Cioq0Lo69s6YAUAgP5/8++8HVZxdu9ib8Cv+vtmzyZs2DZk0iWhtLXX33Xe0BJmdt90/nsPkW+6k8i9/IBgI8PgTT7N+/ftpyeKW3zL7Jq/TmKozrQROF5H+wDZgLDCuHcfmi0ihqtYAw4BVyQ5qcww5Fdo7ZJFubseQzfFJ9RiyyQypGEM+tH2d+zHk3me2+XzxGWW/ITblb56q/lJEZgCrVHVRfAba/wD5wAFgh6qeGT/28BRRAd4CyuMXB48qWQ/ZGGP8JZq60VRVrQQqW2ybnvB4JbGhjNaO/SvQ6sSHo7GCbIzJLNohl7e+EFaQjTGZJfnFOs+ygmyMySzWQzbGGG/Q1M2y+MJZQTbGZJYUXtT7ollBNsZkFhuyMMYYj7CLesYY4xHWQzbGGI+wi3rGGOMRdlHPGGO8Ickthz3NCrIxJrPYGLIxxniEDVkYY4xHWA/ZGGM8wmlId4JjZgXZGJNZbMji6Pz2DRzdn5mX7gjt5rfX2JgOZUMWxhjjET7uISf71mljjPGXaNT9koSIDBeRKhHZKCI/b2X/UBFZLSKNIjKmlf0ni8g2EZnpJrr1kI0xGUVTdFFPRILALOASIAKsFJFFqro+odlHwPXAlKOc5m7gFbfPaT1kY0xm0aj7pW1lwEZV/SD+bdFPAaObPZXqFlV9FzjiZCIyGOgJvOQ2uhVkY0xmSd2QRQmwNWE9Et+WlIgEgIeAn7YnuhVkY0xmaUcPWUTKRWRVwlKecCZp7ewuU/w7UKmqW5O2TGBjyMaYzNKOWRaqOgeYc5TdEaBPwnopsN3lqb8BnC8i/w7kAp1EpF5Vj7gwmMgKsjEms6RuHvJK4HQR6Q9sA8YC41xFUB1/+LGIXA+ck6wYgw1ZGGMyTWOj+6UNqtoITAJeBDYAC1V1nYjMEJFRACJyrohEgCuB2SKy7niiWw/ZGJNZUvhJPVWtBCpbbJue8HglsaGMts7xOPC4m+ezgmyMySw+/qSeFWRjTGaxe1kYY4xHWA/ZGGM8wnrIxhjjEUlmT3iZFWRjTGZRtx+m8x4ryMaYzOLjMWTPfDCkU1kZBfPnU1BRQddxrXwYJjubvOnTKaiooPujjxIoLo5tz8ri5Ntvp/u8eXSfO5fss876/JisLLrddhsFCxZQMH8+OUOHdlj+V19fxcixE7n8qgnMXbDwiP2r3l7DlTdM4qtDv81Lf1/WbF/1jp388JZpXDGunFHjy9lW/XGHZPT7a5zMZZdeyLq1/+C99a/ys5/enLYc7eG3zL7Im8L7IX/RvNFDDgToNnkytVOm4NTU0P23v+Xga6/hfPhhU5MuI0YQra9n9/jx5AwbRrfycvbOmEGXkSMB2DNhAnLKKeTfdx97fvxjUOWka64hWlvL7muvBRHk5JM7JL7jONzz0Cwe+829FBf14PsTJ3PRt77Oaf37NrXp1bOIe+64jcf/+Kcjjp96z4OU/2As3yw7m88+248EWrunyXHy+Wuc/McL8MjDv2T4iKuJRKp5fUUlzz3/Ehs2/CstedzwW2bf5PXxRb129ZBF5FsicquIXJrKENkDB+Js24ZTXQ2NjRxYsoScIUOatckZMoQDL7wAwMFXXqHT4MEAZPXty6HVqwHQ2lqi9fVkhcNArMB8WlERO4EqundvKmM3WbPhfU4t7U2fkl5kZ2dz+cUXsGTZ683alPTqSXhAfwLSvNhu2vwhjuPwzbKzAejatQtdOndOeUa/v8bJlJ37NTZt2sLmzR/R0NDAwoXPMuqKy9KSxS2/ZfZNXsdxv3hMmwVZRN5MePxDYCbQDfiP1r7O5JhDFBYSralpWo/W1BAsLGzWJlhYiHO4jeMQra9H8vJo3LQpVliCQQLFxWSHwwSLipDcXAByJ0yg+5w55N11F4H8/FRFbmZnzS6Kiz7P27OoBztrdrs6dsvWbXTLzWXy1LsZc/3NPDhzLk4H/EXx+2ucTO+SYrZGPr8RV2RbNb17F6cli1t+y+ybvD4eskjWQ85OeFwOXKKqvwAuBca3fgjN7jG6YLvbu9W14OZKqSr7Fy+O/Qo+ezbdJk2iYe3a2DtfMEiwqIiGtWvZU15Ow7p15N5007FlOYao4nLUwXEcVr+zlimTJvLU3EeIbN/Bnyv/ltqAR+Oj1zgZaeUFV49fbfdbZt/k9XFBTjaGHBCRfGKFW1S1BkBVPxWRo072S7zH6McXXpj0TyxaU0MgobcWKCzE2bWrWRsn3qOL1tTEemq5uWhdHQD1s2Y1tcufOZPGSATduxfdv5+Dy2IX0A4sXUr+iBHJohyTnkU92LHz897nxzt3UdijwN2xhT0YGDqNPiW9ABg29Bu8u+49ILW/Cvr9NU5mW6SaPqW9m9ZLS3pR3UEXR1PFb5l9kzeDx5DzgLeAVUB3ESkGEJFcWr+b/jFpqKoiWFoau6qflUXnYcM4uHx5szYHly+n8/DhAORccEHTmCY5ORAfc+00eDA4TtOFqoMrVjTNCOg0eDCNCRewUunLA0N8FNlOZPsOGhoaWPzyK1z0rfPcHXtGiLp99ez5pBaAN996h9P6nZryjH5/jZNZueptBgzoT79+fcjOzuaqq0bz3POuv8osLfyW2S95NaquF69ps4esqv2OsisKfDdlKRyHfQ8/TP4DD0AgwIHFi3G2bOGkG26gsaqKg8uXs7+ykrxp0yioqEDr6tg7YwYAgfx88u+/H1Rxdu1i7733Np123+zZ5E2bhkyaRLS2lrr77ktZ5ERZWUGm/eQmfnTrnTiOw3dHXsqAL/Vl5mPzOXNgiIvOP481G6q4Zerd1O2rZ+lrbzBr7pM8WzGbYDDIlJsncuPkqaAwKDyAMaOGpz6kz1/j5D+ew+Rb7qTyL38gGAjw+BNPs379+2nJ4pbfMvsmrweHItySjh4DcjNk4SXdn5mX7gjttmfMhHRHaLeS5R6bKmU8ofHQtuP+zfuzWZNc15yuN8/sgDmmx84b85CNMSZVfNxDtoJsjMksVpCNMcYjvDgVzyXP3MvCGGNSIoXzkEVkuIhUicjG1j4MJyJDRWS1iDSKyJiE7WeJyAoRWSci74rI991Etx6yMSazpGg6m4gEgVnAJUAEWCkii1R1fUKzj4DrgSktDv8M+IGq/ktEegNviciLqlrb1nNaQTbGZJbU3XqgDNioqh8AiMhTwGigqSCr6pb4vmbdbVV9P+HxdhHZCRQCbRZkG7IwxmQUjUZdL4m3eYgv5QmnKgG2JqxH4tvaRUTKgE7ApmRtrYdsjMks7RiySLzNQytam6PcrvEQEekFLACuU03+mW4ryMaYzJK6e1lEgD4J66WA67ulicjJwF+AO1X19WTtwYYsjDGZJqrul7atBE4Xkf4i0gkYCyxyEyHe/n+A+ar6326jW0E2xmSWRsf90gZVbQQmAS8CG4CFqrpORGaIyCgAETlXRCLAlcBsEVkXP/wqYChwvYi8HV/OauVpmrEhC2NMZknh7TdVtRKobLFtesLjlcSGMloe9yTwZHufzwqyMSazePC2mm5ZQTbGZBS1e1kYY4xHWA/ZGGM8wgpy5vDjzd79eFN9ep+f7gQmU3XAt7Z/UawgG2Myihe/K88tK8jGmMxiBdkYYzzCZlkYY4xHWA/ZGGM8wgqyMcZ4gzo2ZGGMMd5gPWRjjPEGm/ZmjDFeYQXZGGM8wr9DyFaQjTGZRRv9W5GtIBtjMot/67EVZGNMZvHzRT37Tj1jTGaJtmNJQkSGi0iViGwUkZ+3sn+oiKwWkUYRGdNi33Ui8q/4cp2b6NZDNsZklFT1kEUkCMwCLgEiwEoRWaSq6xOafQRcD0xpcWx34D+AcwAF3oof+0lbz2k9ZGNMZkldD7kM2KiqH6jqIeApYHRiA1XdoqrvtnK2y4C/quqeeBH+KzA82RNaQTbGZBRtdL+ISLmIrEpYyhNOVQJsTViPxLe5cUzHeqYgdyoro2D+fAoqKug6btyRDbKzyZs+nYKKCro/+iiB4uLY9qwsTr79drrPm0f3uXPJPuusz4/JyqLbbbdRsGABBfPnkzN06Amf+bBXX1/FyLETufyqCcxdsPCI/aveXsOVN0ziq0O/zUt/X9ZsX/WOnfzwlmlcMa6cUePL2Vb9cYdkbK/LLr2QdWv/wXvrX+VnP7053XFc8VtmP+TVaDsW1Tmqek7CMifhVNLa6V3GOKZjvVGQAwG6TZ5M7e23s/u66+g8bBjBvn2bNekyYgTR+np2jx/Pp888Q7fy2BtZl5EjAdgzYQKfTJlCt5tuAom9Fiddcw3R2lp2X3stu6+7jkPvvHNiZ45zHId7HprFfz10N4sqZlP5t6Vs2vxhsza9ehZxzx23MeKSi444fuo9D3LDuDE894c5PPXYw3TPz0t5xvYKBAI88vAvGXnFNXzlqxfx/e9/hzPOOD3dsdrkt8y+yZu6IYsI0CdhvRTY7jLFMR3b7oIsIvPbe0wy2QMH4mzbhlNdDY2NHFiyhJwhQ5q1yRkyhAMvvADAwVdeodPgwQBk9e3LodWrAdDaWqL19WSFw0CsIH5aURE7gSq6d+8JnfmwNRve59TS3vQp6UV2djaXX3wBS5a93qxNSa+ehAf0JyDN3+g3bf4Qx3H4ZtnZAHTt2oUunTunPGN7lZ37NTZt2sLmzR/R0NDAwoXPMuqKy9Idq01+y+yXvO3pISexEjhdRPqLSCdgLLDIZYwXgUtFJF9E8oFL49va1GZBFpFFLZbngP/78LrLYEkFCguJ1tQ0rUdraggWFjZrEywsxDncxnGI1tcjeXk0btoUK4TBIIHiYrLDYYJFRUhuLgC5EybQfc4c8u66i0B+fqoi+zLzYTtrdlFc9HnWnkU92Fmz29WxW7Zuo1tuLpOn3s2Y62/mwZlzcTzwpZK9S4rZGvm8AxLZVk3v3sVpTJSc3zL7JW+qCrKqNgKTiBXSDcBCVV0nIjNEZBSAiJwrIhHgSmC2iKyLH7sHuJtYUV8JzIhva1OyHnIpUAf8J/BQfNmX8LhViQPlC7a77eG3oC6GalTZv3gxTk0N3WfPptukSTSsXRv71tlgkGBREQ1r17KnvJyGdevIvemmY8uSYZlbiymtjXi1wnEcVr+zlimTJvLU3EeIbN/Bnyv/ltqAx0Ba+QHUzZ9HGvkts1/yqiOul6TnUq1U1ZCqnqaqv4xvm66qi+KPV6pqqaqepKoFqnpmwrHzVHVAfPm9m+zJ5iGfA0wG7gB+qqpvi8h+VX0lyQ8xB5gD8PGFFyb9E4vW1BBI6F0GCgtxdu1q1saJ90CjNTWxnmVuLlpXB0D9rFlN7fJnzqQxEkH37kX37+fgstgFqQNLl5I/YkSyKK75MfNhPYt6sGPn5737j3fuorBHgbtjC3swMHQafUp6ATBs6Dd4d917xGb5pM+2SDV9Sns3rZeW9KLaIxcbj8Zvmf2S18VQhGe12UNW1aiq/hq4AbhDRGbSAR8maaiqIlhaGpuFkJVF52HDOLh8ebM2B5cvp/Pw2DS+nAsuaBqDJScH4mOYnQYPBsfB+TB2gergihVNMxg6DR5M44fNL1ydaJkP+/LAEB9FthPZvoOGhgYWv/wKF33rPHfHnhGibl89ez6pBeDNt97htH6npjxje61c9TYDBvSnX78+ZGdnc9VVo3nu+ZfSHatNfsvsl7waFdeL17gqrqoaAa4UkW8TG8JILcdh38MPk//AAxAIcGDxYpwtWzjphhtorKri4PLl7K+sJG/aNAoqKtC6OvbOmAFAID+f/PvvB1WcXbvYe++9TafdN3s2edOmIZMmEa2tpe6++07szHFZWUGm/eQmfnTrnTiOw3dHXsqAL/Vl5mPzOXNgiIvOP481G6q4Zerd1O2rZ+lrbzBr7pM8WzGbYDDIlJsncuPkqaAwKDyAMaOSznfvcI7jMPmWO6n8yx8IBgI8/sTTrF//frpjtclvmf2S1889ZOnoMSA3Qxbm+HR/Zl66I7Rbl97npzuC8aDGQ9uOu9u67RvDXNeckhVLPNVNtntZGGMyip97yFaQjTEZJepi9oRXWUE2xmQUL16sc8sKsjEmo1hBNsYYj/DgZ1Vcs4JsjMko1kM2xhiPULWCbIwxnuDYLAtjjPEG6yEbY4xH2BiyMcZ4hM2yMMYYj7AesjHGeIQT9cZXhR4LK8jGmIzi5yEL/76VGGNMK6IqrpdkRGS4iFSJyEYR+Xkr+3NE5On4/jdEpF98e7aIPCEia0Rkg4hMdZPdCrIxJqOoiuulLSISBGYBlwODgKtFZFCLZjcCn6jqAODXwOFvlLgSyFHVrwCDgR8dLtZtsYJsjMkoqu6XJMqAjar6gaoeAp4CRrdoMxp4Iv74GeBiiX0brAIniUgW0AU4hItvW+rwMeSS5f/q6KcwPvz2jf3bl6U7QrvsGTMh3RGMS26GIg4TkXKgPGHTnPiXNAOUAFsT9kWAr7c4RVMbVW0Ukb1AAbHiPBqoBroCP1HVPcny2EU9Y0xGac8si3jxnXOU3a1V9pb96qO1KQMcoDeQDywTkb+p6gdt5bEhC2NMRtF2LElEgD4J66XA9qO1iQ9P5AF7gHHAC6raoKo7gdeAc5I9oRVkY0xGSeEsi5XA6SLSX0Q6AWOBRS3aLAKuiz8eAyzR2DdHfwQMk5iTgPOA95I9oQ1ZGGMySqpuLhQfE54EvAgEgXmquk5EZgCrVHUR8DtggYhsJNYzHhs/fBbwe2AtsWGN36vqu8me0wqyMSajpPJLp1W1EqhssW16wuMDxKa4tTyuvrXtyVhBNsZkFG31Ops/WEE2xmSURrsfsjHGeIP1kI0xxiNSOYb8RbOCbIzJKNZDNsYYj7AesjHGeIRjPWRjjPEGH3+DkxVkY0xmiVoP2RhjvMHH3+BkBdkYk1nsop4xxnhEVPw7ZOHL229edumFrFv7D95b/yo/++nN6Y7jit8yezHvq6+vYuTYiVx+1QTmLlh4xP5Vb6/hyhsm8dWh3+alvzf/RpLqHTv54S3TuGJcOaPGl7Ot+uMOydiprIyC+fMpqKig67hxRzbIziZv+nQKKiro/uijBIqLY9uzsjj59tvpPm8e3efOJfussz4/JiuLbrfdRsGCBRTMn0/O0KEnfOa2OO1YvMZ3PeRAIMAjD/+S4SOuJhKp5vUVlTz3/Ets2ODdr4ryW2Yv5nUch3semsVjv7mX4qIefH/iZC761tc5rX/fpja9ehZxzx238fgf/3TE8VPveZDyH4zlm2Vn89ln+5FAB/SiAgG6TZ5M7ZQpODU1dP/tbzn42ms4H37Y1KTLiBFE6+vZPX48OcOG0a28nL0zZtBl5EgA9kyYgJxyCvn33ceeH/8YVDnpmmuI1tay+9prQQQ5+eQTO3MSfp5l4bsectm5X2PTpi1s3vwRDQ0NLFz4LKOuuCzdsdrkt8xezLtmw/ucWtqbPiW9yM7O5vKLL2DJstebtSnp1ZPwgP4EWvzKumnzhziOwzfLzgaga9cudOncOeUZswcOxNm2Dae6GhobObBkCTlDhjRrkzNkCAdeeAGAg6+8QqfBgwHI6tuXQ6tXA6C1tUTr68kKh4FYQfy0oiJ2AlV0794TOnMyUcT14jXHXJBF5IZUBnGrd0kxWyOff4tKZFs1vXsXpyOKa37L7MW8O2t2UVxU2LTes6gHO2t2uzp2y9ZtdMvNZfLUuxlz/c08OHMujpP6X1gDhYVEa2qa1qM1NQQLC5u1CRYW4hxu4zhE6+uRvDwaN22KFcJgkEBxMdnhMMGiIiQ3F4DcCRPoPmcOeXfdRSA//4TOnEwKv8LpC3c8PeRfHG2HiJSLyCoRWRWNfnocT9HquY/Ypi6+zzud/JbZi3lbe3q3124cx2H1O2uZMmkiT819hMj2Hfy58m+pDXg0bl43VfYvXhwbMpg9m26TJtGwdi04DgSDBIuKaFi7lj3l5TSsW0fuTTdZ5jZExf3iNW2OIYvI0b5yRICeRzsu8ZtcszqVpPRf8rZINX1Kezetl5b0orqDLtCkit8yezFvz6Ie7Nj5eU/u4527KOxR4O7Ywh4MDJ1Gn5JeAAwb+g3eXfcekNphmGhNDYGE3mWgsBBn165mbZx4DzRaUxPrWebmonV1ANTPmtXULn/mTBojEXTvXnT/fg4ui12kPLB0KfkjRpzQmZPx87S3ZD3knsAPgCtaWdz9vphiK1e9zYAB/enXrw/Z2dlcddVonnv+pXREcc1vmb2Y98sDQ3wU2U5k+w4aGhpY/PIrXPSt89wde0aIun317PmkFoA333qH0/qdmvKMDVVVBEtLY7MQsrLoPGwYB5cvb9bm4PLldB4+HICcCy5oGoMlJwfi49qdBg8Gx2m6sHZwxYqmGQydBg+mMeGC24mYORlH3C/JiMhwEakSkY0i8vNW9ueIyNPx/W+ISL+Eff8mIitEZJ2IrBGRpBcuks2yeB7IVdW3WwmyNOlP0wEcx2HyLXdS+Zc/EAwEePyJp1m//v10RHHNb5m9mDcrK8i0n9zEj269E8dx+O7ISxnwpb7MfGw+Zw4McdH557FmQxW3TL2bun31LH3tDWbNfZJnK2YTDAaZcvNEbpw8FRQGhQcwZtTw1Id0HPY9/DD5DzwAgQAHFi/G2bKFk264gcaqKg4uX87+ykrypk2joKICratj74wZAATy88m//35Qxdm1i7333tt02n2zZ5M3bRoyaRLR2lrq7rvvxM6cRKp6yCISJPZlpZcAEWCliCxS1fUJzW4EPlHVASIyFrgP+L6IZAFPAteq6jsiUgA0JH3Ojh4bTPWQhckM+7cvS97IQ/aMmZDuCCeEnkuXHvfI7uzSa1zXnB9Fnjzq84nIN4C7VPWy+PpUAFX9VUKbF+NtVsSL8A6gELgcGKeq17Qnu++mvRljTFtU3C+JExDiS3nCqUqArQnrkfg2Wmujqo3AXqAACAEqIi+KyGoR+Zmb7ON9nOkAABPBSURBVL77YIgxxrSlPUMWiRMQWtFa77ll7/tobbKAbwHnAp8BL4vIW6r6clt5rIdsjMkoKfzodATok7BeCmw/Wpv4kEUesCe+/RVV3aWqnwGVwNnJntAKsjEmo6RwHvJK4HQR6S8inYCxwKIWbRYB18UfjwGWaOzC3IvAv4lI13ihvgBYTxI2ZGGMySipmmWhqo0iMolYcQ0C81R1nYjMAFap6iLgd8ACEdlIrGc8Nn7sJyLyn8SKugKVqvqXZM9pBdkYk1FS+cEQVa0kNtyQuG16wuMDwJVHOfZJYlPfXLOCbIzJKH6eZ2sF2RiTUbx4jwq3rCAbYzKKF28875YVZGNMRon6eNDCCrIxJqP4+W5vVpCNMRnFv/1jK8jGmAxjPWRjjPGIRvFvH9kKsjEmo/i3HFtBNsZkGBuyMKad/HbD9+7PzEt3hHbz22ucKjbtzRhjPMK/5dgKsjEmw9iQhTHGeITj4z6yFWRjTEaxHrIxxniEWg/ZGGO8wXrIxhjjEX6e9mZfcmqMySjajiUZERkuIlUislFEft7K/hwReTq+/w0R6ddi/6kiUi8iU9xkt4JsjMkojajrpS0iEgRmAZcDg4CrRWRQi2Y3Ap+o6gDg18B9Lfb/GljsNrsVZGNMRtF2/JdEGbBRVT9Q1UPAU8DoFm1GA0/EHz8DXCwiAiAi3wE+ANa5zW4F2RiTUaLtWESkXERWJSzlCacqAbYmrEfi22itjao2AnuBAhE5Cbgd+EV7sttFPWNMRmnPtDdVnQPMOcru1r4uteXJj9bmF8CvVbU+3mF2xQqyMSajpHDaWwTok7BeCmw/SpuIiGQBecAe4OvAGBG5HzgFiIrIAVWd2dYTWkE2xmQUR1M27W0lcLqI9Ae2AWOBcS3aLAKuA1YAY4AlqqrA+YcbiMhdQH2yYgxWkI0xGSZV85BVtVFEJgEvAkFgnqquE5EZwCpVXQT8DlggIhuJ9YzHHs9zWkE2xmSUVH50WlUrgcoW26YnPD4AXJnkHHe5fT4ryMaYjOLnj077ctrbZZdeyLq1/+C99a/ys5/enO44rvgts1fydioro2D+fAoqKug6ruXwHZCdTd706RRUVND90UcJFBfHtmdlcfLtt9N93jy6z51L9llnfX5MVhbdbruNggULKJg/n5yhQzss/6uvr2Lk2IlcftUE5i5YeMT+VW+v4cobJvHVod/mpb8va7avesdOfnjLNK4YV86o8eVsq/64QzL6/TVuKYq6XrzGdz3kQCDAIw//kuEjriYSqeb1FZU89/xLbNjwr3RHOyq/ZfZM3kCAbpMnUztlCk5NDd1/+1sOvvYazocfNjXpMmIE0fp6do8fT86wYXQrL2fvjBl0GTkSgD0TJiCnnEL+ffex58c/BlVOuuYaorW17L72WhBBTj65Q+I7jsM9D83isd/cS3FRD74/cTIXfevrnNa/b1ObXj2LuOeO23j8j3864vip9zxI+Q/G8s2ys/nss/1IwP30Kdd8/hq3xs93e0vaQxaRgSJyu4g8IiIPxx+f8UWEa03ZuV9j06YtbN78EQ0NDSxc+CyjrrgsXXFc8Vtmr+TNHjgQZ9s2nOpqaGzkwJIl5AwZ0qxNzpAhHHjhBQAOvvIKnQYPBiCrb18OrV4NgNbWEq2vJyscBmIF5tOKitgJVNG9ezsk/5oN73NqaW/6lPQiOzubyy++gCXLXm/WpqRXT8ID+hNoMVd10+YPcRyHb5adDUDXrl3o0rlzyjP6/TVujaPqevGaNguyiNxO7OOCArxJbBqIAH9s7UYbX4TeJcVsjXw+FTCyrZrevYvTEcU1v2X2St5AYSHRmpqm9WhNDcHCwmZtgoWFOIfbOA7R+nokL4/GTZtihSUYJFBcTHY4TLCoCMnNBSB3wgS6z5lD3l13EcjP75D8O2t2UVz0ed6eRT3YWbPb1bFbtm6jW24uk6fezZjrb+bBmXNxHCflGf3+GrfGz0MWyXrINwLnqur/VtUn48v/JvYZ7xuPdlDixxGj0U9TmZfWPvWiHnynS+S3zJ7O6yaHKvsXL479Cj57Nt0mTaJh7VpwHAgGCRYV0bB2LXvKy2lYt47cm276wqK6/dCW4zisfmctUyZN5Km5jxDZvoM/V/4ttQGPxkevcWva89Fpr0k2hhwFegMfttjeizZ+nsSPI2Z1Kknpv+RtkWr6lPZuWi8t6UV1B13sSBW/ZfZK3mhNDYGE3lqgsBBn165mbZx4jy5aUxPrqeXmonV1ANTPmtXULn/mTBojEXTvXnT/fg4ui11AO7B0KfkjRnRI/p5FPdix8/Pe58c7d1HYo8DdsYU9GBg6jT4lvQAYNvQbvLvuPSC1Q0d+f41bk8ljyLcAL4vIYhGZE19eAF4GJnd8vCOtXPU2Awb0p1+/PmRnZ3PVVaN57vmX0hHFNb9l9krehqoqgqWlsav6WVl0HjaMg8uXN2tzcPlyOg8fDkDOBRc0jWmSkwPxMddOgweD4zRdqDq4YkXTjIBOgwfT+GHL/kZqfHlgiI8i24ls30FDQwOLX36Fi751nrtjzwhRt6+ePZ/UAvDmW+9wWr9TU57R769xa/w8ZNFmD1lVXxCRELEhihJi48cRYKWqpn5AywXHcZh8y51U/uUPBAMBHn/iadavfz8dUVzzW2bP5HUc9j38MPkPPACBAAcWL8bZsoWTbriBxqoqDi5fzv7KSvKmTaOgogKtq2PvjBkABPLzyb//flDF2bWLvffe23TafbNnkzdtGjJpEtHaWurua3kL29TIygoy7Sc38aNb78RxHL478lIGfKkvMx+bz5kDQ1x0/nms2VDFLVPvpm5fPUtfe4NZc5/k2YrZBINBptw8kRsnTwWFQeEBjBk1PPUhff4at8Yzw2vHQDo6fKqHLExm2PbN09MdoV26PzMv3RHabc+YCemO0G49ly497rl9l/YZ7rrmvLT1hQ6YS3jsfDcP2Rhj2uLFoQi3rCAbYzKKn4csrCAbYzKK9ZCNMcYj/DztzQqyMSajePEj0W5ZQTbGZBQbsjDGGI+wgmyMMR7h51kWvrxBvTHGHE0qPzotIsNFpEpENrZ2h0sRyRGRp+P73xCRfvHtl4jIWyKyJv7/YW6yW0E2xmQUbcd/bRGRIDALuBwYBFwtIoNaNLsR+ERVBwC/Bg5/RnwXcIWqfoXYt1IvcJPdCrIxJqM4GnW9JFEGbFTVD1T1ELF7w49u0WY08ET88TPAxSIiqvpPVT18U/F1QGcRyUn2hFaQjTEZRVVdL0mUAFsT1iPxba22UdVGYC/Q8h6r3wP+qaoHkz2hXdQzxmSU9syyEJFyoDxh05z4/dwhdnfLllqevM02InImsWGMS93ksYJsjMko7fmkXuKXabQiAvRJWC8Fth+lTUREsoA8YA+AiJQC/wP8QFU3ucljQxbGmIwSVXW9JLESOF1E+otIJ2AssKhFm0XELtoBjAGWqKqKyCnAX4Cpqvqa2+xWkI0xGSVVsyziY8KTgBeBDcBCVV0nIjNEZFS82e+AAhHZCNwKHJ4aNwkYAPwvEXk7vhQly25DFsaYjOJi9oRrqloJVLbYNj3h8QHgylaOuwe4p73PZwXZGBf8+O0bfvyWk1RwMRThWVaQjTEZxW6/aYwxHmE9ZGOM8QjrIRtjjEc46qQ7wjGzgmyMySh+vv2mFWRjTEaxG9QbY4xHWA/ZGGM8wmZZGGOMR9gsC2OM8YhUfnT6i2YF2RiTUWwM2RhjPMLGkI0xxiOsh2yMMR5h85CNMcYjrIdsjDEeYbMsjDHGI/x8Uc+X36l32aUXsm7tP3hv/av87Kc3pzuOK37L7JW8ncrKKJg/n4KKCrqOG3dkg+xs8qZPp6Cigu6PPkqguDi2PSuLk2+/ne7z5tF97lyyzzrr82Oysuh2220ULFhAwfz55AwdesLmbenV11cxcuxELr9qAnMXLDxi/6q313DlDZP46tBv89LflzXbV71jJz+8ZRpXjCtn1PhytlV/3GE526Kqrhev8V1BDgQCPPLwLxl5xTV85asX8f3vf4czzjg93bHa5LfMnskbCNBt8mRqb7+d3dddR+dhwwj27dusSZcRI4jW17N7/Hg+feYZupWXx7aPHAnAngkT+GTKFLrddBOIAHDSNdcQra1l97XXsvu66zj0zjsnZt4WHMfhnodm8V8P3c2iitlU/m0pmzZ/2KxNr55F3HPHbYy45KIjjp96z4PcMG4Mz/1hDk899jDd8/M6JGcyqfqSUwARGS4iVSKyUUR+3sr+HBF5Or7/DRHpl7Bvanx7lYhc5iZ70oIsIgNF5GIRyW0Z1M0TpFrZuV9j06YtbN78EQ0NDSxc+CyjrnD1s6aN3zJ7JW/2wIE427bhVFdDYyMHliwhZ8iQZm1yhgzhwAsvAHDwlVfoNHgwAFl9+3Jo9WoAtLaWaH09WeEwECuKn1ZUxE6giu7de0LmbWnNhvc5tbQ3fUp6kZ2dzeUXX8CSZa83a1PSqyfhAf0JxN8sDtu0+UMcx+GbZWcD0LVrF7p07twhOZNJVQ9ZRILALOByYBBwtYgMatHsRuATVR0A/Bq4L37sIGAscCYwHHg0fr42tVmQReT/BZ4F/h9grYiMTth9b7KTd4TeJcVsjWxvWo9sq6Z37+J0RHHNb5m9kjdQWEi0pqZpPVpTQ7CwsFmbYGEhzuE2jkO0vh7Jy6Nx06ZYMQwGCRQXkx0OEywqQnJj/YrcCRPoPmcOeXfdRSA//4TM29LOml0UF32et2dRD3bW7HZ17Jat2+iWm8vkqXcz5vqbeXDmXBwnPTeKj6q6XpIoAzaq6geqegh4Chjdos1o4In442eAi0VE4tufUtWDqroZ2Bg/X9uSvHusAXLjj/sBq4DJ8fV/tnFcebztKqC8Pe9YLpYrVXXu4edR1WtV9f9L8XOkerlSVecmvBZez+yV17gpRxuv2zpVLT28/sknn+xU1QJVzVLVX6vq26r6rKpWqupoVe2hMd+LH3Orqi5IV15V3TRo0KCfpClvsyUUCl0ZCoXmJqxfGwqFWv1zP/XUU5eHQqExCW3HhEKhvaFQ6EuhUCgrFAr9KRQK3djBfz+Oe2lRq5rVK2AMMDdh/VpgZovj1wKlCeubgB7ATOCahO2/A8Yky5NsyCKoqvXxwr0FuBC4XET+E5CjHaSqc1T1nPgyJ8lztFcE6BN/XA6UAtuP3twTDmcuj697PbNXXuPEHBwlR2KbLFXNB/YAjcBPgLOI9VZOAf4F7AY+A/4nfsx/A2enKy+Qt2HDhvFpytuSm/wANDQ0hFo59p9VVVUfVFVVNQJ/puNypkyLWtWyXrVW41p2q4/Wxs2xR0hWkHeISNPl3nhxHknsHeAryU7eQVYCpwP9c3JyhNg4zaI0ZXFrJXB6OBzuBHTC+5m98ho35eDor9si4Lr44zErVqzYR+wvflfgpPj2S4gVvPXxfc8R61wAXBzfnpa8wJJ4DyodeVs6/Pe0f/zvanv+3FcC+eFw+PCYxzA6LucXpV1vsCKSBeQR6xC4fnNrJkl3vhQoPsq+IWn8VWOEqr7/0UcfHVDVO9L9a4/bzJs3bz6gqpt8ktkrr/EIVX2/xes2Q1VHxR93VtX/VtWNqvrmGWec8W58ez9VrVLVDar6N1Xtm3DOvqr6D1V9V1VfVtVT05VXVb8ErEpj3mZLKBQaEQqF3g+FQptCodAd8W0zQqHQqPjjc0OhUGTAgAFOKBTaHQqF1iUce0koFHo3FAqtCYVCj4dCoU5p/Htz3Aux32A+4PM32HeAM1u0uRn4bfzxWGBh/PGZ8fY58eM/IDbi0PZzpvuHPs4XLNXj05bX55n9ltcye3sBRgDvExsbviO+bQYwKv64M7FhpI3Am8CXEo69I35cFXC5m+eT+IHGGGPSzHcfDDHGmExlBdkYYzzCdwVZRDqLyJsi8o6IrBORX6Q7kxsiskVE1ojI2yKyKt153BCRoIj8U0SeT3eWZERknojsFJG16c7SHsk+muslItJHRP4uIhvi//YmpztTpvHdGHL8UzAnqWq9iGQDrxL7sMrrSQ5NKxHZApyjqrvSncUtEbkVOAc4WVVHpjtPW0RkKFAPzFfVL6c7jxvxj9K+T2yaW4TY1LGrVdWT08VEpBfQS1VXi0g34C3gO17N60e+6yFrTH18NTu++OtdxQdEpBT4NjA33VncUNV/EJv/6SduPprrGaparaqr44/3ARuAkvSmyiy+K8jQ9Kv028BO4K+q+ka6M7mgwEsi8paIlCdtnX6/AX4G+Pdu395XAmxNWI/gkwIXv6vZ1wA//NvzDV8WZFV1VPUsYh9cKRMRP/yKOkRVzyZ256ib479ie5KIjAR2qupb6c6S4Y7p47XpFr/z45+AW1S1Lt15MokvC/JhqloLLCV2eztPU9Xt8f/vJHZfguR3fkqfIcCo+Lj3U8AwEXkyvZEy0rF9vDaN4tdt/gRUqOr/SXeeTOO7giwihSJySvxxF+D/At5Lb6q2ichJ8YsgiMhJwKXE7hLlSao6VWN3JOtH7OOgS1T1mjTHykQrgdNFpL+IeP4eJ/EL6r8DNqjqf6Y7TybyXUEGegF/F5F3if2F/quqen1aVk/gVRF5h9jHK/+iqi+kOVNGEZE/AiuAsIhEROTGdGdKRlUbgUnAi8QukC1U1XXpTdWmIcRuQTksPn3zbREZke5QmcR3096MMSZT+bGHbIwxGckKsjHGeIQVZGOM8QgryMYY4xFWkI0xxiOsIBtjjEdYQTbGGI/4/wHQRjO4apT1zgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Number of blocks (clusters in a matrix)\n",
    "num_blocks = 2\n",
    "\n",
    "# The size of each block (cluster) in a modeled matrix\n",
    "block_size = 3\n",
    "\n",
    "# Correlation inside a block\n",
    "block_corr = 0.6\n",
    "\n",
    "# Correlation between the clusters\n",
    "std = 0.4\n",
    "\n",
    "# For the same result output\n",
    "np.random.seed(3)\n",
    "\n",
    "# Finding the random vector of means and covariance matrix\n",
    "mu_vec, cov_matrix = nco.form_true_matrix(num_blocks, block_size, block_corr, std)\n",
    "\n",
    "#Plotting the heatmap\n",
    "sns.heatmap(cov_matrix, annot=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The result is a covariance matrix with exactly 2 blocks (clusters) that have 3 elements each.\n",
    "\n",
    "As the function shuffles the resulting covariance matrix, in other outputs of the function the block structure may not be seen so clearly as in this example. \n",
    "\n",
    "This covariance matrix and vector of means can be used to run the MCOS algorithm."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Number of observations to use to obtain empirical covariance matrix and the means vector\n",
    "num_obs = 100\n",
    "\n",
    "# Number of simulations to do in the MCOS\n",
    "num_sims = 50\n",
    "\n",
    "# The bandwidth of the KDE kernel.\n",
    "kde_bwidth = 0.25\n",
    "\n",
    "# Flag if the goal is the minimum variance\n",
    "min_var_portf = True\n",
    "\n",
    "# Flag if we want to use the Ledoit-Wolf shrinkage procedure\n",
    "lw_shrinkage = False\n",
    "\n",
    "# For the same result output\n",
    "np.random.seed(1)\n",
    "\n",
    "# Finding the optimal weights for minimum variance\n",
    "w_cvo, w_nco = nco.allocate_mcos(mu_vec, cov_matrix, num_obs, num_sims, kde_bwidth, min_var_portf, lw_shrinkage)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The result is a DataFrame of weigths of NCO and CVO for every simulation. We can make boxplots to examine the differences in the results."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAskAAAF1CAYAAAAa1Xd+AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3df5xddX3g/9fbhAAr1KLYyC8BK7oTYouSYndNNWNQcGHBXbEygKXrVJb9mtQt7leg42JFp1/S7mJbpIvo8JVWGbTa2ihUanWGmoe/wIraMEUjRYhoqeKvCAKJ7/3jnMGbkzszZ5J7c++5eT0fj/vIPed8zrnve+fmPe/5nM/5nMhMJEmSJP3ME3odgCRJktRvLJIlSZKkCotkSZIkqcIiWZIkSaqwSJYkSZIqLJIlSZKkCotkdUREvCci3tbB470tIr4TEd/u1DGbKiI2R8Samm3viYiTuxySpIYzZ3dPRPxaRNxVs+2aiNja7Zi0eyySB0xZJD0cEdsi4nsRcVNEHNXruFpFREbEM+fZfhTwBmBFZj6tQ68ZEfHbEfGPEfHjiNgaEX8REc+JiEsj4u/b7HNoRDwaESvL5SMj4n0R8d3yGJ+PiNM7Ed98MvP4zJze0+OYjKX+Y86e9zW/EhFPaFn3toh4T8vysoj4vYj4WpmT74mI6yLimJY2p5e5+sdl7n5fRBzZiRjnkpmfysxnd+JYnf5jRotjkTyY/mNmHgQcBvwLcFWP41mso4HvZuYDi90xIpbOsemPgdcDvw08GXgW8GHgNODPgX8fEcdW9jkb+Epm/mNEPBnYBDwKHA8cCrwduCEizlpsnJLUwpzd3uEUeXguHwTOAM4BngT8MvAFYG157LOAGyjy/6EUufsRYFNEHLLYWLXvsUgeYJn5E4oksmJ2XUQ8KSL+LCL+NSK+ERFvmv1LPSL+T0R8sKXthoj4RNkLu6bsff3d8pTaPRFx7lyvHRGvjYgtEfFgRGyMiMPL9bM9tl8qe05eVdnvZODjwOHl9veU688ohx18PyKmI2KoZZ97IuLiiPgy8ONq0o2I44DXASOZ+cnMfCQzH8rM92XmFZm5Ffgk8OrK2/gN4Pry+e8A24DRzPx2Zj6cmZPAOPC/IyLafAbXR8QbyudHlD0j/0+5/Mzys4ly+fSIuKN8f5+OiF+qvL+Ty+cHlsf9XkTMRMQb2/QOnxARX46IH0TE+yPigIh4IvA3LZ/rtog4PCJOiojbI+KHEfEvEXFl2x+opK4zZ+/iD4C3tNtevu5LgDMz87bM3J6ZP8jMqzNzosyt/xt4W5nrH87MbwO/RZHLf6fNMQ+Iolf/0HL5TRGxPSJ+rlx+W0T8Ufl8/4j4XxFxb5k7r4mIA8ttO521i4jnRcQXI+JHUZzBfH9Ueocj4g0R8UBEfCsi/ku57gLgXOCN5Wf7kXL9xRHxzfJ4d0XE2jk+P+2pzPQxQA/gHuDk8vm/oSjy/qxl+58Bfw0cDBwDfJWi8Jtt/1XgN4FfA74DHFluWwNsB64E9gdeBPwYeHa5/T0UyQjgxeW+zyvbXgX8fUsMCTxznvewBtjasvys8rVeAuwHvBHYAixrec93AEcBB7Y53oXANxb43M4Fvtay/GyKXuOnlsufBd7SZr9jy/fz7DbbXgN8pHx+DvB14P0t2/66fP484AHg+cAS4PzyPe3f5md6BXArcAhwJPDlymd1D/B5ih6YJwMzwIXtPtdy3WeAV5fPDwJ+tdffYR8+9qUH5uxdcnbLax5H0TP8W+W6twHvKZ9fAdw6T0z/tjzGsW22vQX4zBz7/T3wivL531Lk7Ze1bPtP5fM/AjaWefZg4CPA/1f9PIBlwDcozmTuB/xnit8tb2tpux24vNz+H4CHgEOqP6dy+dnAfcDh5fIxwC/2+ns8qA97kgfThyPi+8APKZLUHwJExBLgVcClmfmjzLyH4i/tVwNk5kPAeRRJ9b3A+ix6WVv9zyx6Ym8FbgJ+vc3rnwtcl5n/kJmPAJcC/y5axokt0quAmzLz45n5GPC/gAOBf9/S5k8y877MfLjN/k8BvrXAa/wVsDwiZo/5G8DfZOa/lsuHznGMb7Vsr7oV+LWy1+eFFL0iLyi3vajcDvBa4J2Z+bnM3JGZ11OcEvzVNsf8deD3M/N75c/mT9q0+ZPMvD8zH6RI3Ce0f8sAPAY8MyIOzcxtmfnZedpK6g5zdnsJ/E/gsojYv7Jtobw+m5PnytvtcjYUeflFZe/1L1Hk2BdFxAHArwCfKnupXwv8TmY+mJk/An6f9kNDfhVYSvF+H8vMv6ToyGj1GHB5uf1mip7uucY076D4Q2ZFROyXmfdk5tfnaKs9ZJE8mF6emT9P8R9pHXBrRDyNIinM/lU76xvAEbMLmfl54G4ggA9Ujvu9zPxxZd/D27z+4a2vkZnbgO+2vs4iVY/3U4q/pFuPd988+3+XYqzfnMpfNn8B/EaZAM/lZ0MtoOhlaXeMw1q2V4/5dYpkdwJFL89Hgfsj4tnsXCQfDbyhPC35/fKX5VHM/dm2vtd277v16vKHKHqI5zJK0evzTxFxW+yFCxEl7cKcPYeyaLwXuKCyaaG8PpuT58rbu+Ts0q0UvbvPA75CMZTkRRTF7pbM/A7wVIpe/C+05OyPleurDge+mVl0+5aq7/27mbm9ZXnOvJ2ZW4D/Dvwe8EBE3Dg7NEadZ5E8wMpeyb+k+MtzNUVSeIyiKJv1dOCbswsR8TqKRH0/xSmyVoeUY1tb972/zUvf3/oa5T5PaX2dRaoeLyiKyNbjZXWnFp8AjoyIVQu8zvUUvSwvoTh99tGWbX8HvCJarrQu/TpFwvvqHMe8FTiL4jTjN8vl36AYLnFH2eY+YDwzf77l8W+yGPNc9S2KYRazFnMV/C6fUWZ+LTNHgF8ANgAfrPyMJe0l5uw5vQkYoyhMZ/0dcFLMPVPFXcBW4JWtK8sc/gqK3wvtfJqiF/c/UQznuJPiczuNn3VsfAd4GDi+JWc/KYuLL6u+BRxRfgaz9jRv35CZqyk+46TI3eoCi+QBVl68cSZFQTaTmTsoehrGI+LgiDgauIjiNB0R8SyKMV/nUZzOe2NEVE/VvyWKaXd+DTidove16gbgv0TECeUpst8HPleeKoTi6u1nLOKtfAA4LSLWRsR+FFMNPUKRzBaUmV8D/hSYLC+oWFZeoHF2RFzS0vRTwPeBa4EbM/PRlm1vB34OmIiIp5X7j1Ak7v+30kvQ6laKnqHZi1+mgfXApvLnAfAu4MKIeH75M3tiRJwWEQfP8VlcGhGHRMQR5bHr+hfgKRHxpNkVEXFeRDy17On5frl6R9u9JXWVObu9LKbA/ArF9Rqz6/6Oopf3ryLixIhYWn5GF0bEa8qc/D+AN0XEOVFc9Pw04N0Uufztc7zWQxTjoF/Hz4riTwP/dXa5zJfvAt4eEb8Aj1+cfUqbQ36GIqeuK2M8EzhpEW9/p88+Ip4dES8uf04/oSjWzdldYpE8mD4SEdsoxreNA+dn5uZy23qKCyruppjS7AbgunL81XuBDZn5pbKw/F3gz1vGgn0b+B5FL8H7KC4I+6fqi2fmJyjGkX2I4q/oX2TnsVq/B1xfnqZqNz6uery7KH4JXEXxF/x/pJgy6dF5d9zZbwPvAK6mKAa/TtFT8JGW10mKi2SOLv9tjeG7FD07BwB3Upzqu4jiorf3z/O6t1L0Ss8WyZsoekMen5c5M2+nGN/2DorPdwvFhTjtXE7RO/LPFD0pH6T45bOg8mc1CdxdfvaHA6cCm8vvyx8DZ2dxhb2kvcecvbA3UVwk1+os4Gbg/cAPgH8EVlHkRsrc/GqKmSy+Q5G7DwReUOb0udxKcRHd51uWW/M4wMUUufqzEfHD8jV3GUdcvuf/TDG07fsUn8tHqZm3gQmK8cffj4gPU5w1uKJ8P9+mOAv4uzWPpUWKuTvApJ+J4o5v783Mrk7CrsWJiP9GUdi+qNexSOof5uz+FRGfA67JzP+/17FofvYkSw0SEYdFxAsi4gnlBYBvoJiZQ5LUhyLiReUwvaURcT7FrBkf63VcWlitIjkiTo1iwuotlTGc1XZnRXHDhFUt6y4t97trjvE6kupbBrwT+BHFDVD+mmK8tfQ4c7bUV54NfIliSMgbgLMyc6FpSdUHFhxuEcU8jV+luOJ/K3AbxZ3L7qy0O5hiDsZlwLrMvD0iVlCMgTyJYhqUvwOe1XLBkiSpg8zZktQZdXqST6KYG/DucgD6jcCZbdq9leJmCa0X/ZxJMUvAI5n5zxSD3BdzVackaXHM2ZLUAXWK5CPYeeLrrVQmGI+I5wJHZWbrvLK19pUkdZQ5W5I6YGmNNtFm3eNjNKKYmPvttJ+yat59W45xAeXddA488MATjzpqMfNs99ZPf/pTnvCEZlz/2KRYoVnxNilWaFa8TYoV4Ktf/ep3MrPdnbf2lq7n7PI4jczbTfs+NSneJsUKzYq3SbFCs+KdL2fXKZK3svPdYY5k5zv2HAysBKajuKHM04CNEXFGjX0ByMxrKW7gwKpVq/L222+vEVZ/mJ6eZs2aNb0Oo5YmxQrNirdJsUKz4m1SrAAR8Y2FW3VV13M2NDdvN+371KR4mxQrNCveJsUKzYp3vpxdp8y/DTguIo6NiGUUE4xvnN2YmT/IzEMz85jMPAb4LHBGeYOEjcDZEbF/RBwLHMfPJueWJHWeOVuSOmDBnuTM3B4R64BbgCXAdZm5OSIuB27PzI3z7Ls5Ij5AcZeb7cDrvEpakrrHnC1JnVFnuAWZeTPFrR9b1102R9s1leVxittsSpL2AnO2JO25ZoyqliRJkvYii2RJkiSpwiJZkiRJqrBIliRJkioskiVJkqQKi2RJkiSpwiJZkiRJqrBIliRJkioskiVJkqQKi2RJkiSpwiJZkiRJqrBIliRJkioskiVJkqQKi2RJkiSpwiJZkiRJqrBIliRJkioskiVJkqQKi2RJkiSpwiJZkiRJqrBIliRJkioskiVJkqQKi2RJkiSpwiJZkiRJqrBIliRJkioskiVJkqQKi2RJkiSpwiJZkiRJqrBI3k2Tk5OsXLmStWvXsnLlSiYnJ3sdkiRJkjpkaa8DaKLJyUnGxsaYmJhgx44dLFmyhNHRUQBGRkZ6HJ0kSZL2lD3Ju2F8fJyJiQmGh4dZunQpw8PDTExMMD4+3uvQJEmS1AEWybthZmaG1atX77Ru9erVzMzM9CgiSZIkdZJF8m4YGhpi06ZNO63btGkTQ0NDPYpIkiRJnVSrSI6IUyPirojYEhGXtNl+YUR8JSLuiIhNEbGiXH9MRDxcrr8jIq7p9BvohbGxMUZHR5mammL79u1MTU0xOjrK2NhYr0OTJHO2JHXAghfuRcQS4GrgJcBW4LaI2JiZd7Y0uyEzrynbnwFcCZxabvt6Zp7Q2bB7a/bivPXr1zMzM8PQ0BDj4+NetCep58zZktQZdWa3OAnYkpl3A0TEjcCZwOMJNzN/2NL+iUB2Msh+NDIywsjICNPT06xZs6bX4UjSLHO2JHVAZM6fGyPiLODUzPytcvnVwPMzc12l3euAi4BlwIsz82sRcQywGfgq8EPgTZn5qTavcQFwAcDy5ctPvPHGG/fwbe0927Zt46CDDup1GLU0KVZoVrxNihWaFW+TYgUYHh7+Qmau6tXr742cXe7fyLzdtO9Tk+JtUqzQrHibFCs0K955c3ZmzvsAXgm8u2X51cBV87Q/B7i+fL4/8JTy+YnAfcDPzfd6J554YjbJ1NRUr0OorUmxZjYr3qbEesMNN+Txxx+fT3jCE/L444/PG264odchLagpn+0s4PZcIK9287G3c3Y2LG837fvUpHibFGtms+JtUqyZzYp3vpxdZ7jFVuColuUjgfvnaX8j8H/KAvwR4JHy+Rci4uvAs4Dba7yupA7yJjj7DHO2JHVAndktbgOOi4hjI2IZcDawsbVBRBzXsnga8LVy/VPLi0iIiGcAxwF3dyJwSYvjTXD2GeZsSeqABXuSM3N7RKwDbgGWANdl5uaIuJyii3ojsC4iTgYeA74HnF/u/kLg8ojYDuwALszMB7vxRiTNz5vg7BvM2ZLUGXWGW5CZNwM3V9Zd1vL89XPs9yHgQ3sSoKTOmL0JzvDw8OPrvAnOYDJnS9Ke84570j7Cm+BIklRfrZ5kSc3nTXAkSarPIlnah3gTHEmS6nG4hSRJklRhkSxJkiRVWCRLkiRJFRbJkiRJUoVFsiRJklRhkSxJkiRVWCRLkiRJFRbJkiRJUoU3E5EWKSIW1T4zuxSJJEnqFnuSpUXKzF0eR1/80bbrLZAlSWomi2RJkiSpwuEWkvqGQ1kkSf3CnmRJfcOhLJKkfmFPck32cEmSJO077Emuaa6erLl6uSRJktRcFsmSJElShcMt1HMOZZEkSf3GnmT1nENZJElSv7EnWRpg9tJLkrR77EmWBpi99JIk7R6LZEmSJKnCIlmSJEmqsEiWJEmSKiySJUmSpAqLZEmSJKnCIlmSJEmqsEiWJA289evXc8ABBzA8PMwBBxzA+vXrex3SvCYnJ1m5ciVr165l5cqVTE5O9jokaZ/jzUQkSQNt/fr1XHPNNWzYsIEVK1Zw5513cvHFFwNw1VVX9Ti6XU1OTjI2NsbExAQ7duxgyZIljI6OAjAyMtLj6KR9hz3JkqSB9q53vYsNGzZw0UUXccABB3DRRRexYcMG3vWud/U6tLbGx8eZmJhgeHiYpUuXMjw8zMTEBOPj470OTdqn1CqSI+LUiLgrIrZExCVttl8YEV+JiDsiYlNErGjZdmm5310RcUong5ck7cqcvbNHHnmECy+8cKd1F154IY888kiPIprfzMwMq1ev3mnd6tWrmZmZ6VFE0r5pwSI5IpYAVwMvA1YAI60JtXRDZj4nM08A/gC4stx3BXA2cDxwKvCn5fEkSV1gzt7V/vvvzzXXXLPTumuuuYb999+/RxHNb2hoiE2bNu20btOmTQwNDfUoImnfVKcn+SRgS2benZmPAjcCZ7Y2yMwftiw+Ecjy+ZnAjZn5SGb+M7ClPJ4kqTvM2RWvfe1rufjii7nyyiv5yU9+wpVXXsnFF1/Ma1/72l6H1tbY2Bijo6NMTU2xfft2pqamGB0dZWxsrNehSfuUOhfuHQHc17K8FXh+tVFEvA64CFgGvLhl389W9j2izb4XABcALF++nOnp6Rph9Y9+j/cTn/gE733ve7n33nt5+tOfznnnncfatWt7HVYt/f7ZtmpSrNCseJsUax/oes4u929M3n7FK17B1q1bueSSS3jsscfYb7/9OP3003nFK17Rl3EfdthhnHvuubzmNa/ZKW8fdthhfRnvrG3btvV1fFVNirdJsULz4p1LnSI52qzLXVZkXg1cHRHnAG8Czl/EvtcC1wKsWrUq16xZUyOsPvGxm+jneCcnJ3nf+97Hddddt9NV0itWrOj/q6T7/LPdSZNihWbF26RY+0PXc3a5f6Py9mx809PTjfg+rVmzhre+9a2NiRea89nOalK8TYoVmhfvXOoMt9gKHNWyfCRw/zztbwRevpv7qsPGx8c555xzWL9+Paeccgrr16/nnHPO8SppaXCZsyX1xKDN712nJ/k24LiIOBb4JsVFHee0NoiI4zLza+XiacDs843ADRFxJXA4cBzw+U4ErnruvPNOHnrooV3m27znnnt6HZqk7jBnS9rrBnF+7wV7kjNzO7AOuAWYAT6QmZsj4vKIOKNsti4iNkfEHRRj3M4v990MfAC4E/gY8LrM3NGF96E5LFu2jHXr1u003+a6detYtmxZr0OT1AXmbEm9MIjze9e6415m3gzcXFl3Wcvz18+z7zjQ3E+o4R599FGuuuoqnvvc57Jjxw6mpqa46qqrePTRR3sdmqQuMWdL2tsGcX5vb0s94FasWMHLX/5y1q9fz8zMDENDQ5x77rl8+MMf7nVokiRpQMzO7z08PPz4uqbP722RPODGxsbajhFq8ukPSZLUX2bn956tN2bn925yvWGRPOBmB8u39iSPj483dhC9JEnqP4NYb1gk7wNGRkYYGRkZmHkLJUlS/xm0eqPOPMmSJEnSPsUiWZIkSaqwSJYkSZIqLJIlSZL60KDd5rlpvHBPkiQtSkQsqn1mdimSwTWIt3luGnuSJUnSomTmLo+jL/5o2/UWyLtnEG/z3DQWyZIkSX1mEG/z3DQWyZIkSX1m9jbPrZp+m+emsUiWJEnqM7O3eZ6ammL79u2P3+Z5bGys16HtM7xwT5Ikqc8M4m2em8YiWZIkqQ8N2m2em8bhFpIkSX3IeZJ7y55kSZKkPuM8yb1nT7IkSVKfcZ7k3rNIliRJ6jPOk9x7DrcYUN4yVJKk5pqdJ3l4ePjxdc6TvHfZkzygvGWoJDVHRLR9DA8Pt12vwec8yb1nT7IkST02V2fFMZfcxD1XnLaXo1EvVf8IevGLX7zT8jnnnMM555zz+LIdXd1jT7IkSVKf8Exw/7AnWZI0kLw2Q9KesCdZkjSQ7JGTtCcskiVJkqQKi2RJkiSpwiJZkiRJqrBIliRJkioskiVJkqQKp4CTJEkDy6kAtbvsSZYkSQNrrin/5poOUJpVq0iOiFMj4q6I2BIRl7TZflFE3BkRX46IT0TE0S3bdkTEHeVjYyeDlyTtypwtSXtuweEWEbEEuBp4CbAVuC0iNmbmnS3NvgisysyHIuK/AX8AvKrc9nBmntDhuCVJbZizJakz6vQknwRsycy7M/NR4EbgzNYGmTmVmQ+Vi58FjuxsmJKkmszZktQBdS7cOwK4r2V5K/D8edqPAn/TsnxARNwObAeuyMwPV3eIiAuACwCWL1/O9PR0jbD6R5PibVKs0Kx4mxQrNCveJsXaB7qes6HZebtJsUKz4m1SrNCseJsU67Zt2xoV71zqFMntLgttO7I9Is4DVgEvaln99My8PyKeAXwyIr6SmV/f6WCZ1wLXAqxatSrXrFlTJ/b+8LGbaEy8TYoVmhVvk2KFZsXbpFj7Q9dzNjQ4bzft+9SkeJsUKzQr3j6NddBnDqkz3GIrcFTL8pHA/dVGEXEyMAackZmPzK7PzPvLf+8GpoHn7kG8kqT5mbMl7RWDPnNInSL5NuC4iDg2IpYBZwM7XfEcEc8F3kmRbB9oWX9IROxfPj8UeAHQevGIJKmzzNmS1AELDrfIzO0RsQ64BVgCXJeZmyPicuD2zNwI/CFwEPAXZdf7vZl5BjAEvDMifkpRkF9RucJaktRB5mxJ6oxad9zLzJuBmyvrLmt5fvIc+30aeM6eBChJWhxztiTtOe+4J0mSJFVYJEuSJEkVtYZbSJ3yy2/5W37w8GO12x9zyU212j3pwP340ptfurthSZIk7cQiWXvVDx5+jHuuOK1W2+np6drzQtYtpiVJkupwuIUkSZJUYZEsSZIkVVgkS5IkSRUWyZIkSVKFRbIkSZJUYZEsSZIkVVgkS5IkSRUWyZIkSVKFRbIkSZJUYZEsSZIkVVgkS5IkSRUWyZIkSVKFRbIkSZJUYZEsSZIkVVgkS5IkSRUWyZIkSVKFRbIkSZJUsbTXAUiSpP70y2/5W37w8GO12x9zyU212j3pwP340ptfurthSXuFRbIkSWrrBw8/xj1XnFar7fT0NGvWrKnVtm4xLfWSwy0kSZKkCnuSpQHhaVFJkjrHIlmax2IKz14XnZ4WlSSpcyySpXnULTwtOiVJGiyOSZYkSZIqLJIlSZKkCotkSZIkqcIxyZL2OmfikCT1u1pFckScCvwxsAR4d2ZeUdl+EfBbwHbgX4HXZOY3ym3nA28qm74tM6/vUOySGsqZOLrLnC1Je27B4RYRsQS4GngZsAIYiYgVlWZfBFZl5i8BHwT+oNz3ycCbgecDJwFvjohDOhe+JKmVOVuSOqNOT/JJwJbMvBsgIm4EzgTunG2QmVMt7T8LnFc+PwX4eGY+WO77ceBUYHLPQ+8eTwVLarB9LmdLUjfUKZKPAO5rWd5K0cswl1Hgb+bZ94jFBNgLngqW1GD7XM6WpG6oUyRHm3XZtmHEecAq4EWL2TciLgAuAFi+fDnT09M1wuquujFs27ZtUfH2+r31+vUXE0O/fLZ1jtukWKE/4m1SrA3T9Zxd7tt3ebuuJsUKvY93kP+v9kMMdTUpVmhevG1l5rwP4N8Bt7QsXwpc2qbdycAM8Ast60aAd7YsvxMYme/1TjzxxOy1oy/+aO22U1NTXTluN/T69RcbQz98tnWP26RYM3sfb5NiXSzg9lwgr3bzsbdzdvZJ3q6rH74ji9HreAf5/2o/xFBXk2LNbFa88+XsOvMk3wYcFxHHRsQy4GxgY2uDiHhumUzPyMwHWjbdArw0Ig4pL/54ablOktQd5mxJ6oAFh1tk5vaIWEeRKJcA12Xm5oi4nKL63gj8IXAQ8BcRAXBvZp6RmQ9GxFspkjbA5VleECJJ6jxztiR1Rq15kjPzZuDmyrrLWp6fPM++1wHX7W6AkqTFMWdL0p7zttSSJElShbelliQ1mnPbq4n83vY/i2TtVQcPXcJzrr+k/g41b4h78BBAvbmtJQ0W57ZXE/m97X8WydqrfjRzhUlBkiT1PcckS5IkSRX2JEuSJGle3RhD3e/jpy2SB8BivrgO/JckSYvVjTHU/T5U0iJ5ANT94jrGV5I0yPbF3k51j0WyJEl7kVN/dc++2Nup7rFIliRpL3LqL6kZnN1CkiRJqrBIliRJkioskiVJkqQKi2RJkiSpwgv3JElSWwcPXcJzrr+k/g7X1z0uQL2LF6VesUiWJElt/WjmCmfi0D7L4RaSJElShUWyJEmSVOFwC0mSpL3M8d79zyJZkiRpL3O8d/9zuIUkSZJUYZEsSZIkVTjcQhoQjm+TJKlzLJKlAdGk8W0W9JKkfmeRLGmva1JBL0naN1kkS5IazTMTkrrBIlmS1GiemdCsbvzB5B9LhX3xs7VIluaxqKRg75Qk9VQ3/mDyj6XCvvjZWiRL86ibFOydkiRpsFgkt+H4NkmSpH2bRXIbjm+TJEnat1kkS5KkOS2qg+dj9do+6cD9djMaae+pVSRHxKnAHwNLgHdn5hWV7S8E/gj4JeDszPxgy7YdwFfKxXsz84xOBC5Jas+crU6pe1YVimJ6MU17hsEAABBISURBVO2lfrdgkRwRS4CrgZcAW4HbImJjZt7Z0uxe4DeB/9HmEA9n5gkdiFVzcAYGSbPM2ZLUGXV6kk8CtmTm3QARcSNwJvB4ws3Me8ptP+1CjFqAMzBIamHOlqQOqFMkHwHc17K8FXj+Il7jgIi4HdgOXJGZH642iIgLgAsAli9fzvT09CIO3x11Y9i2bdui4u3We6tz3CbFCs2Kt0mxQn/E26RYG6brORv6L2837fvUtHib8vqLiWExn63fg8Udtx8+247IzHkfwCspxrTNLr8auGqOtu8BzqqsO7z89xnAPcAvzvd6J554Yvba0Rd/tHbbqamprhx3Meoet0mxZjYr3ibFmtn7eJsU62IBt+cCebWbj72ds7MP8nbTvk9Ni7cpr7/YGOp+tn4PFn/cXn+2izFfzn5CjTp6K3BUy/KRwP2LKMLvL/+9G5gGnlt3X0nSopmzJakD6hTJtwHHRcSxEbEMOBvYWOfgEXFIROxfPj8UeAEt4+IkSR1nzpakDliwSM7M7cA64BZgBvhAZm6OiMsj4gyAiPiViNhKcZrvnRGxudx9CLg9Ir4ETFGMbzPhSlKXmLMlqTNqzZOcmTcDN1fWXdby/DaKU3rV/T4NPGcPY5QkLYI5W5L2XJ3hFpIkSdI+xdtSS5K0Fy3qBlDgTaAGmLf87m8WyZIk7UV1bwAF3gRqkHnL7/7ncAtJkiSpwp5kSVLjedpaUqdZJEuSGs3T1pK6weEWkiRJUoVFsiRJklRhkSxJkiRVWCRLkiRJFRbJkiRJUoWzWwyI2tMf9cHUR07VJElS83T693e//+62SB4Adacz6oepj5yqqbv8A0SS1A374u9vi2RpQOyLCUySqva13k51j0WyJEkaCHYWqJO8cE+SJEmqsEiWJEmSKiySJUmSpAqLZEmSJKnCIlmSJEmqcHYLSZL2Muc0l/qfRbIkSXuR05RJzWCRLC2gSbf8liRJnWGRLM2jSbf8bhpPN0uS+plFsqS9ztPNkqR+5+wWkiRJUoVFsiRJklThcIs5OF5SkiRp32WR3IbjJSVJkvZtDreQJEmSKiySJUmSpAqLZEmSJKmiVpEcEadGxF0RsSUiLmmz/YUR8Q8RsT0izqpsOz8ivlY+zu9U4JKk9szZkrTnFiySI2IJcDXwMmAFMBIRKyrN7gV+E7ihsu+TgTcDzwdOAt4cEYfsediSpHbM2ZLUGXV6kk8CtmTm3Zn5KHAjcGZrg8y8JzO/DPy0su8pwMcz88HM/B7wceDUDsQtSWrPnC1JHVBnCrgjgPtalrdS9DLU0W7fI6qNIuIC4AKA5cuXMz09XfPw/aFJ8TYpVmhWvE2KFZoVb5Ni7QNdz9nQjLw9PDzcdn1saN9+amqqi9Hsvn78bOfSpFihWfE2KVZoXrzt1CmSo826rHn8Wvtm5rXAtQCrVq3KNWvW1Dx8H/jYTTQm3ibFCs2Kt0mxQrPibVKs/aHrORuakbczdw19enq6Wd+nJn3/mxQrNCveJsUKzYt3DnWGW2wFjmpZPhK4v+bx92RfSdLimbMlqQPqFMm3AcdFxLERsQw4G9hY8/i3AC+NiEPKiz9eWq6TJHWHOVuSOmDBIjkztwPrKBLlDPCBzNwcEZdHxBkAEfErEbEVeCXwzojYXO77IPBWiqR9G3B5uU6S1AXmbEnqjDpjksnMm4GbK+sua3l+G8VpuXb7XgdctwcxSpIWwZwtSXvOO+5JkiRJFbV6kiVJktR9Ee0mmZl76sJ2s7ioM+xJliRJixIRuzy+seH0tuvnKvrUXmbu8piammq73gK5uyySJUnSoljIaV9gkSxJkiRVWCRLkiRJFV64p56bb7xauwsVPHUnSZK6zZ5k9dxcY9jmGt8mSZLUbRbJkiRJUoVFsiRJklRhkSxJkiRVWCRLkiRJFRbJkiRJUoVTwEmSpIHlNKPaXfYkS5KkgeU0o9pdFsmSJElShUWyJEmSVGGRLEmSJFVYJEuSJEkVFsmSJElShUWyJEmSVGGRLEmSJFVYJEuSJEkVFsmSJElShUWyJEmSVGGRLEmSJFVYJEuSJEkVFsmSJElShUWyJGngTU5OsnLlStauXcvKlSuZnJzsdUg7iYi2j29sOL3teqkfDPr3dmmvA5AkqZsmJycZGxtjYmKCHTt2sGTJEkZHRwEYGRnpcXSFzGy7fnp6mjVr1uzdYKSaBv17a0+yJGmgjY+PMzExwfDwMEuXLmV4eJiJiQnGx8d7HZqkPmaRLEkaaDMzM6xevXqndatXr2ZmZqZHEUlqglpFckScGhF3RcSWiLikzfb9I+L95fbPRcQx5fpjIuLhiLijfFzT2fAlSVXm7J0NDQ2xadOmndZt2rSJoaGhHkUk1dPvY+kH3YJjkiNiCXA18BJgK3BbRGzMzDtbmo0C38vMZ0bE2cAG4FXltq9n5gkdjlsLmGuAfGxo336ucUWSmsWcvauxsTFGR0cfH5M8NTXF6Oiowy3U15owln7Q1elJPgnYkpl3Z+ajwI3AmZU2ZwLXl88/CKyNJl7GOEAyc5fH1NRU2/UWyNJAMWdXjIyMMD4+zvr16znllFNYv3494+PjFhrqa46l7706s1scAdzXsrwVeP5cbTJze0T8AHhKue3YiPgi8EPgTZn5qeoLRMQFwAUAy5cvZ3p6ejHvoeeaEu+2bdsaEys0L94mxQrNirdJsfaBrudsaF7ePuyww3jHO97Btm3bOOigg4BmfK+alAebFCv0f7wzMzPs2LGD6enpx2PdsWMHMzMzfR039P9nW1edIrld70K163GuNt8Cnp6Z342IE4EPR8TxmfnDnRpmXgtcC7Bq1aps1LQhH7upMdOcNG1KlkbF26DvAdCseJsUa3/oes6G5ubtRuUVmhVvk2KF/o93aGiIJUuWsGbNmsdjnZqaYmhoqK/jhv7/bOuqM9xiK3BUy/KRwP1ztYmIpcCTgAcz85HM/C5AZn4B+DrwrD0NWpI0J3O2NABmx9JPTU2xffv2x8fSj42N9Tq0fUadnuTbgOMi4ljgm8DZwDmVNhuB84HPAGcBn8zMjIinUiTeHRHxDOA44O6ORS9JqjJnSwNgdsz8+vXrmZmZYWhoyLH0e9mCRXI5Xm0dcAuwBLguMzdHxOXA7Zm5EZgA/jwitgAPUiRlgBcCl0fEdmAHcGFmPtiNNyJJMmdLg2RkZISRkZGBGb7QNLVuS52ZNwM3V9Zd1vL8J8Ar2+z3IeBDexijJGkRzNmStOe8454kSZJUYZEsSZIkVVgkS5IkSRUWyZIkSVKFRbIkSZJUYZEsSZIkVVgkS5IkSRUWyZIkSVKFRbIkSZJUYZEsSZIkVVgkS5IkSRUWyZIkSVKFRbIkSZJUYZEsSZIkVSztdQBNERFzb9uw67rM7GI0kiRJ6iZ7kmvKzLaPqamptuslSZLUXBbJkiRJUoVFsiRJklRhkSxJkiRVWCRLkiRJFRbJ+4DJyUlWrlzJ2rVrWblyJZOTk70OSZIkqa85BdyAm5ycZGxsjImJCXbs2MGSJUsYHR0FYGRkpMfRSZIk9Sd7kgfc+Pg4ExMTDA8Ps3TpUoaHh5mYmGB8fLzXoUmSJPUti+QBNzMzw+rVq3dat3r1amZmZnoUkSRJUv+zSB5wQ0NDbNq0aad1mzZtYmhoqEcRSZIk9T+L5AE3NjbG6OgoU1NTbN++nampKUZHRxkbG+t1aJIkSX3LC/cG3OzFeevXr2dmZoahoSHGx8e9aE+SJGkeFsn7gJGREUZGRpienmbNmjW9DkeSJKnvOdxCkiRJqrBIliRJkioskiVJkqSKWkVyRJwaEXdFxJaIuKTN9v0j4v3l9s9FxDEt2y4t198VEad0LnRJC4mIto9vbDi97XoNBnO2JO25BYvkiFgCXA28DFgBjETEikqzUeB7mflM4O3AhnLfFcDZwPHAqcCflseTtBdkZtvH1NRU2/VqPnO2JHVGnZ7kk4AtmXl3Zj4K3AicWWlzJnB9+fyDwNoouqXOBG7MzEcy85+BLeXxJGkXi+n1tud7TuZsSeqAOkXyEcB9Lctby3Vt22TmduAHwFNq7itJQPue77l6ve35npM5W5I6oM48ye26a6q/neZqU2dfIuIC4IJycVtE3FUjrn5xKPCdXgdRU5NihWbFe2hsaEys0LDPlubECnB0j1+/6zkbGp23m/Z9alK8TYoVmhVvk2KFZsU7Z86uUyRvBY5qWT4SuH+ONlsjYinwJODBmvuSmdcC19aIpe9ExO2ZuarXcdTRpFihWfE2KVZoVrxNirVPdD1nQ3PzdtO+T02Kt0mxQrPibVKs0Lx451JnuMVtwHERcWxELKO4qGNjpc1G4Pzy+VnAJ7M4F7oROLu8kvpY4Djg850JXZLUhjlbkjpgwZ7kzNweEeuAW4AlwHWZuTkiLgduz8yNwATw5xGxhaI34uxy380R8QHgTmA78LrM3NGl9yJJ+zxztiR1Rp3hFmTmzcDNlXWXtTz/CfDKOfYdB8b3IMZ+16TTjU2KFZoVb5NihWbF26RY+4I5e15N+z41Kd4mxQrNirdJsULz4m0rvEJckiRJ2pm3pZYkSZIqLJJ300K3fe0nEXFdRDwQEf/Y61gWEhFHRcRURMxExOaIeH2vY5pPRBwQEZ+PiC+V8b6l1zEtJCKWRMQXI+KjvY5lIRFxT0R8JSLuiIjbex2Pms283R1Nytvm7O4atJztcIvdUN6m9avASyimTLoNGMnMO3sa2Bwi4oXANuDPMnNlr+OZT0QcBhyWmf8QEQcDXwBe3sefbQBPzMxtEbEfsAl4fWZ+tsehzSkiLgJWAT+Xmaf3Op75RMQ9wKrMbMp8m+pT5u3uaVLeNmd316DlbHuSd0+d2772jcz8e4or2PteZn4rM/+hfP4jYIY+vuNXFraVi/uVj779yzMijgROA97d61ikvcy83SVNytvmbC2GRfLu8date0FEHAM8F/hcbyOZX3kq7A7gAeDjmdnP8f4R8Ebgp70OpKYE/jYivlDe4U3aXebtvaAJeduc3VUDlbMtkndP7Vu3avdExEHAh4D/npk/7HU888nMHZl5AsXdyU6KiL48NRoRpwMPZOYXeh3LIrwgM58HvAx4XXkKWtod5u0ua0reNmd31UDlbIvk3VP71q1avHKc2IeA92XmX/Y6nroy8/vANHBqj0OZywuAM8oxYzcCL46I9/Y2pPll5v3lvw8Af0VxylzaHebtLmpi3jZnd96g5WyL5N1T57av2g3lRRUTwExmXtnreBYSEU+NiJ8vnx8InAz8U2+jai8zL83MIzPzGIrv7Ccz87wehzWniHhieREQEfFE4KVA31/pr75l3u6SJuVtc3b3DGLOtkjeDZm5HZi97esM8IHM3NzbqOYWEZPAZ4BnR8TWiBjtdUzzeAHwaoq/mO8oH/+h10HN4zBgKiK+TPFL+OOZ2ffT9DTEcmBTRHwJ+DxwU2Z+rMcxqaHM213VpLxtzu6egcvZTgEnSZIkVdiTLEmSJFVYJEuSJEkVFsmSJElShUWyJEmSVGGRLEmSJFVYJEuSJEkVFsmSJElShUWyJEmSVPF/ATu9strxwIXzAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 864x432 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Creating subplot figure with having two side by side plots\n",
    "fig, axes = plt.subplots(nrows=1, ncols=2, figsize=(12,6))\n",
    "\n",
    "# Plot the boxplot of CVO weigths\n",
    "w_cvo.boxplot(ax = axes[0]) \n",
    "axes[0].set_title('Boxplot for CVO weights')\n",
    "\n",
    "# Setting the same grid size\n",
    "axes[0].set_xlim(0, 7)\n",
    "axes[0].set_ylim(0, 0.4)\n",
    "\n",
    "# Plot the boxplot of NCO weigths\n",
    "w_nco.boxplot(ax = axes[1])\n",
    "axes[1].set_title('Boxplot for NCO weights')\n",
    "\n",
    "# Setting the same grid size\n",
    "axes[1].set_xlim(0, 7)\n",
    "axes[1].set_ylim(0, 0.4)\n",
    "\n",
    "# Plot the result\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The weigths produced by the NCO are grouped at higher density near the mean values (have a lower variance of weights).\n",
    "\n",
    "This would mean that the NCO has produced more robust estimates of the optimal weights allocation."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can check it by calculating the mean of the standard deviation between the weights of the method and the true optimal allocation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The mean of the standard deviation with the true weight for NCO is:\n",
      "0.04621200914715077\n",
      "The mean of the standard deviation with the true weight for CVO is:\n",
      "0.06902384964077275\n"
     ]
    }
   ],
   "source": [
    "# Finding the errors in estimations\n",
    "err_cvo, err_nco = nco.estim_errors_mcos(w_cvo, w_nco, mu_vec, cov_matrix, min_var_portf)\n",
    "\n",
    "# Outputting the result\n",
    "print('The mean of the standard deviation with the true weight for NCO is:')\n",
    "print(err_nco)\n",
    "print('The mean of the standard deviation with the true weight for CVO is:')\n",
    "print(err_cvo)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Indeed, the error of the NCO is much lower than the error from simple CVO."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Conclusion"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This notebook describes the Nested Clustered Optimization (NCO) algorithm, the Monte Carlo Optimization Selection (MCOS) algorithm alongside the De-noising algorithm and other helping functions. Also, it shows how these can be used on some real examples.\n",
    "\n",
    "The algorithms and the descriptions were originally presented by _Marcos Lopez de Prado_ in the paper __A Robust Estimator of the Efficient Frontier__  [available here](https://papers.ssrn.com/abstract_id=3469961).\n",
    "\n",
    "Key takeaways from the notebook:\n",
    "- Convex optimization solutions in the problems of portfolio optimization tend to be unstable.\n",
    "- The instability comes from two sources\n",
    "  - Noise in the input variables\n",
    "  - The signal structure that magnifies the estimation errors in the input variables\n",
    "- The De-noising algorithm calculates the eigenvalues of the correlation matrix and eliminates the ones that are higher than the theoretically estimated ones, as they are caused by noise.\n",
    "- The Convex Optimization Solution (CVO) can either find optimal allocations for maximum Sharpe ratio or minimum variance.\n",
    "- The Nested Clustered Optimization (NCO) breaks the correlation matrix into clusters and applies CVO individually to each cluster as well as to the correlation matrix where each element represents a cluster.\n",
    "- The Monte Carlo Optimization Selection (MCOS) allows running optimization methods (in our case CVO and NCO) on multiple sets of empirical covariation matrices and empirical vectors of means. It can show which method is more robust.\n",
    "- The Sample Data Generating function can be used to get a covariance matrix that has distinct clusters with correlations inside a cluster and different correlations between clusters."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "mlfinlab",
   "language": "python",
   "name": "mlfinlab"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
